The unrestricted Skyrme-tensor time-dependent Hartree-Fock 
and its application to the nuclear response from spherical to triaxial nuclei 
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The nuclear time-dependent Hartree-Fock model formulated in the three-dimensional space, based 
on the full Skyrme energy density functional and complemented with the tensor force, is presented 
for the first time. Full self-consistency is achieved by the model. The application to the isovector 
giant dipole resonance is discussed in the linear limit, ranging from spherical nuclei ( 16 O, 120 Sn) to 
systems displaying axial or triaxial deformation ( 24 Mg, 28 Si, 178 Os, 190 W , 238 U). 
Particular attention is paid to the spin-dependent terms from the central sector of the functional, 
recently included together with the tensor. They turn out to be capable of producing a qualitative 
change on the strength distribution in this channel. The effect on the deformation properties is also 
discussed. The quantitative effects on the linear response are small and, overall, the giant dipole 
energy remains unaffected. 

Calculations are compared to predictions from the (quasi)-particle random phase approximation 
and experimental data where available, finding good agreement. 

PACS numbers: 21.60.Jz, 24.30. Cz, 21.60.Ev, 21.30.Fe 
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I. INTRODUCTION 

This work presents the full implementation of the 
Skyrme energy density functional (EDF), including 
the tensor terms, in time-dependent Hartree-Fock 
(TDHF) 0. The TDHF equations have been formu- 
lated on a three dimensional (3D) Cartesian grid with 
no spatial symmetry restrictions, and include all time- 
odd terms. It is consequently possible to study different 
phenomena in both intrinsically spherical and deformed 
systems including triaxiality, within the same formalism 
and without neglecting terms in the p-h channel of the 
nuclear (and Coulomb) interaction. They can range from 
nuclear structure features like giant resonances (includ- 
ing non-linearities and decay by particle emission) and 
large amplitude dynamics like collisions of nuclei (includ- 
ing fusion and fission processes). In this paper we focus 
on the first application, while another one about nucleus- 
nucleus collisions [2| will follow as an extension of recent 
related work [|-Q7 

The Skyrme-tensor force, included in the original de- 
scription of the Skyrme interaction [5[, had mostly been 
omitted from Hartree-Fock calculations after initial ex- 
plorations 0. It has been the object of considerable in- 
terest in the last years, when particular effort has been 
spent to attest the real limits of the mean-field approach. 
For a review of the "new generation" of tensor studies 
we refer to [3] , the authors of which extensively analysed 
the performance of several Skyrme-tensor functionals on 
ground-state properties of spherical systems. In such 
cases, only the so-called time-even terms of the EDF, 



built on densities which maintain the same sign under 
the time-reversal operation, contribute. Their coupling 
constants are typically fixed through a selection of exper- 
imental and empirical data. 

A situation similar to the spherical case, although not 
identical, occurs in axially deformed nuclei. The full 
spatial rotational invariance is lost in the intrinsic sys- 
tem of reference, giving rise to extra contributions to the 
one-body potential with respect to the spherical case [1] ; 
however, the parity and the total angular momentum 
projection along the symmetry axis are still good quan- 
tum numbers, the time-reversal invariance is preserved 
and all the time-odd terms of the functional still van- 
ish in the ground-state. This fact is maintained when 
pairing effects are introduced, as long as the particle- 
particle scattering involves pairs in time-reversal states. 
In [9|], the performance of time-even tensor terms fixed 
on spherical systems Q was studied against deformed 
ground-state properties. Less analyses are available for 
triaxial ground-states, where the spatial symmetry prop- 
erties are further broken, but the time-reversal invariance 
is still preserved (cf. [l(|- II | and references therein). 
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Increasing the complexity of the problem, a different 
situation occurs in odd-A and odd-odd nuclei and, in 
general, in presence of cranking, for which cases the time- 
reversal invariance is broken and the Kramer's degener- 
acy is removed [Hj]-[l4]]- In such situations, in fact, time- 
odd terms become active. 

The effect from the time-odd terms that originate from 
the Skyrme central and spin-orbit force have been studied 
through their influence on rotational bands of supcrdc- 
formed nuclei [l5j . More recently, the full Skyrme time- 
odd sector including the tensor terms has been taken into 
account to investigate the moments of inertia associated 
to superdeformed rotational bands [l6j . 
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Regarding vibrational states, on the side of the (non- 
rclativistic) self-consistent random phase approximation 
(RPA), which does not model the pairing correlations 
and which corresponds to the small amplitude (linear) 
limit of TDHF, progress has been obtained by includ- 
ing the residual Skyrme-tensor interaction in formula- 
tions on a spherical p-h basis. As known, the time- 
odd terms play an essential role in explaining the nu- 
clear response in some cases. Calculations including the 
tensor have been published in both some "neutral" [I?} 
and charge-exchange channels [l8j . Similarly to previous 
work, the relative contribution of the residual two-body 
versus the underlying mean-field tensor effects was dis- 
cussed by these authors, but no further considerations 
on the various functional's terms and their interplay was 
made. Skyrme-tensor RPA calculations in infinite mat- 
ter have been performed as well. 

No self-consistent implementation of the Skyrme func- 
tional complemented by the tensor interaction in a de- 
formed (Q)RPA is currently available. Our task with the 
more general TDHF approach is also aimed at filling this 
gap, allowing more applications on the same footing. In 
this TDHF model all the Skyrme-tensor EDF terms have 
been derived and implemented; the specific role played by 
some of them, where the focus is on spin-dependent terms 
that were previously neglected, will be outlined and is- 
sues related to self-consistency discussed. Although the 
comparison between TDHF and RPA approaches is not 
the aim of this paper and more can be found elsewhere, 
in Section [II A II we will recall some (standard) terminol- 
ogy and elements useful to follow the discussion. 
Previous applications of TDHF to zero temperature gi- 
ant resonances, including time-odd terms from modern 
Skyrme EDFs and exploiting computational power which 
was not available in the past, exist. In [20j , the lin- 
ear response in some light spherical and deformed nuclei 
was studied by comparing different approaches (TDHF 
with absorbing boundary conditions and continuum RPA 
formulated with the Green's functions formalism); the 
terms depending on the square of the momentum den- 
sity (j 2 ), which affects the effective mass, and on the 
square of the spin density (S 2 ) were discussed, while 
other spin-dependent terms were omitted from those cal- 
culations. The authors of [2l| pointed out the relevance 
of the Skyrme time-odd spin-orbit in suppressing spuri- 
ous spin excitations in free translational motion, associ- 
ated to a non physical energy dissipation in heavy-ion 
collisions. This was explained as a consequence of the 
Galilean invariance restoration when the spin-orbit sec- 
tor of the functional is fully included. 

After the formalism is introduced (Sec.|TT|, we present 
(Sec. IIIip the results about the isovector dipole response 
(IVD) for some representative cases. It is well known that 
most of the IVD transition strength is concentrated in 
the giant dipole resonance (IVGDR or simply GDR) , the 
first collective nuclear excitation that was discovered [22T ]- 
[24j . Although the major effect from the inclusion of 
the spin- dependent terms is not expected on this reso- 



nance, we focus on it according to the tradition of the 
TDHF calculations, as this is the simplest, experimen- 
tally well known, non spherical excitation that can be 
reproduced. Many theoretical simulations are available, 
ranging from more phenomenologic approaches (where 
the important dependence on external input can produce 
accurate calculations, but drastically limits the predic- 
tive power), down to microscopic descriptions with no 
ad-hoc adjustable parameters, to which group our model 
belongs. 

In Sec. IIVI conclusions will be summarised. The deriva- 
tion of the tensor contribution to the Skyrme EDF is 
provided in App. [2] as an extension of [2fj, where the 
functional form of the energy density from the central and 
spin-orbit part of the force was written without assuming 
a time-reveral invariant system. Finally, the expressions 
and possible implementations of the densities and cur- 
rents of the functional have been added in App. [Bj 



II. FORMALISM 

The first part of this Section recalls the basic features 
of the TDHF theory and the numerics adopted in this 
work, while the second one focuses on the full Skyrme- 
tensor energy density functional, implemented in TDHF 
with no symmetry restrictions. Some of the performed 
verification and validation tests will be discussed. 



A. Time-dependent Hartree-Fock model 

As mentioned in the Introduction, the theoretical 
framework chosen to implement the unrestricted Skyrme 
energy density functional, to which the following Sec- 
tion is dedicated, is the 3D Cartesian time-dependent 
Hartree-Fock. In this theory 26], with semi-classic limit 
given by the Vlasov equation, well known in plasma 
physics and astrophysics, the interaction among parti- 
cles is modelled in terms of the interaction with the 
one-body potential, which, after (or under) the action 
of an external field, changes in time through the de- 
pendence on the density itself. The master expression 
of the theory derives from the Von Neumann equation, 
the quantum-mechanical version of the classical Liouville 
equation which provides the time evolution of the A-body 
density matrix. Through the constraint p\ = pA for 
the full density matrix, the system can be represented 
by a pure state, solution of the associated Schroedinger 
equation. In TDHF, the previous relation is assumed to 
hold for the one-body density matrix p\ (hereafter p), so 
that the A-particle wave-function can be represented by 
a Slater determinant at any time (the backward relation 
holds too; details can be found in [13|). As a matter of 
fact, the Von Neumann equation, which gives rise to the 
BBGKY hierarchy for the reduced 1-body, 2-body, etc. 
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densities, now simplifies to a one-body problem 



ihStp 



h,p 



(1) 



which can be equivalently represented by the (non-linear) 
set of one-body equations 



ih5 t 4> {i) (x,t) = h[p{t)}^\x,t), 



(2) 



one for each spinor i/j^(x,t) = J2u=2m 4>^ l \ x ^ ljJ ^)i{ ljJ ) 
identified by the index i, with initial conditions given 
further on (the dependence of the one-body Hamilto- 
nian on the particle density can be generalised). A time- 
dependent external field can be inserted in the previous 
equation ([2]). In this work, it will be assumed to simply 
act as a delta function in time, used to prepare the ini- 
tial conditions for the time evolution. The external field 
acting on the system can be defined, in particular, as a 
one-body operator able to induce an isovector response 
(protons and neutrons move in opposition of phase) with 
no charge- mixing, that is in the form 



(3) 



where a+ and a s are creation and annihilation operators 
in the HF basis and F — Df(xi,uji)T z (i). This ex- 
ternal field transfers energy (regulated by the strength 
D) and a selection of angular momentum to the nucleus, 
causing a displacement of the proton and neutron cen- 
tres of mass, compatible deformation and/or spin fluc- 
tuations; in such a way, the system is set in oscillation 
with all the possible frequencies. In other words, after 
the external action, the A-particle wave-function is not 
an eigenstate of the static ( "unperturbed" ) Hamiltonian 
any more and starts to evolve in time. In a microscopic 
picture, the external kick induces p-h transitions that 
are allowed by the selection rules, so that, at time i=0, 



each Hartree Fock eigenstate 
wave-packet 



is transformed into a 



^ (x, t = 0) = e iP 4>P (x) = J2 aatr ( x )> ( 4 ) 

i 

where the index I spans the full HF spectrum. The label 
(i) at the l.h.s. of the previous equation simply enu- 
merates the wave-functions obtained by boosting the HF 
solution having a set of quantum numbers i. Using Q 
as initial condition, where only the choice of the external 
operator F is arbitrary and (partly) defines the problem 
under study, the equations <j2j) are solved and a new set of 
quasiparticles {^^(x,t)} is given at any time. Clearly, 
under the action of the one-body potential, mixings other 
than those initially generated by F are produced. 

According to long-standing prescriptions [28] , in prac- 
tical implementations the time is discretised and the evo- 
lution operator is implemented in a Taylor expansion. 
The time step size and the order at which the expansion 
is arrested are two mutually dependent parameters cho- 
sen in order to ensure the total energy and norms are 



(0) , 



conserved in the "dynamics" within acceptable accuracy. 
The system response is analysed by following the time- 
evolution of the expectation value for the one-body op- 
erator of interest 



(d)(t) = (*(x,t)\d\*(x,t)) 



(5) 



(referred to the ground-state value). One can replace 
the Slater determinant $(X, t) (X — x%, x%, ■ ■ ■) by the 
simple product of the involved TDHF wave- functions and 
recast the previous equation as 



(6)(t) = ^2(^(x,t)\0(x)\^(x,t)), 



(6) 



which can be conveniently expressed in terms of properly 
defined one-body densities. For example, the well known 
(effective) isovector dipole response along the A direction 
is based on the expectation values 



(Of- 01 >(*) 



D 



A {Z Pn (x, t) - N Pp {x, t)) dx, (7) 



depending on the time-dependent isovector density 
Piq{x, t) — p n (x, t) — p p (x, t). It is obtained by choosing 
f(x,co) to be 



fiK{x,uj) = ^iic(»)|cc|I CT 



(8) 



where I a is the identity operator in the spin space 
(kept implicit in what follows) and {D^iif }k=-i,o,1 = 

1 \ ItF TxT ( * s ^ ne se t of / = 1 real spherical har- 

monies written in Cartesian coordinates (D = D\ 

V Y 4-7T ' 

e=l). The subscript in Eq. [7] denotes the J T ST quan- 
tum numbers. The signal obtained on top of a spherical 
ground-state is clearly invariant under any rotation in 
space of the boost direction. The centre of mass correc- 
tion in the definition of the operator is included in order 
to remove the translational zero mode from the strength 
function (the boost needs proper weigths as well). 

Under the hypothesis O = F, it is possible to show 
that the transition strength distribution associated to O 



S 6 {E) = Y,MO\-)\ 2 5{E-E u ), 



(9) 



where v labels the excited states in the small amplitude 
limit, can be obtained through 



S 6 (E) = — 



7T F[Dg(t)} ' 



(10) 



where J- denotes the Fourier transform operation and 
g(t) is the generalization of the (factorized) time profile 
of the external operator. The proof of the equivalence 
between Eqs. © and (10]) can be found in [lj|-|3(J (the 
rearrangement arising from density-dependent Hamilto- 
nian terms is usually considered negligible). 
The presence of the factor D in the denominator of (TIT 
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guarantees that the response is kept constant for differ- 
ent boost strength in the small amplitude (RPA) regime. 
This rescaling can be used as an indicator to test the 
validity of the assumption of linearity. In fact, when 
the boost amplitude is increased above a certain thresh- 
old, the channel for non-harmonic effects, which, under a 
one-body view, accounts for 3p-lh (3h-lp) vertices [3lj ]. 
is opened. 

Choices different from O = F allow one to select the sig- 
nal of interest, like a specific multipole or spin (isospin) 
component, among the various response channels which 
can be opened by a given external operator and mixed up 
during the dynamic evolution. For example, the 1~ states 
are a mixture of both non spin-flip (5=0, 5=1) and spin- 
flip (5=1) components, that is, in the microscopic picture 
of the RPA, nucleons can also reverse the spin when un- 
dergoing particle-hole transitions. A large fraction of the 
isovector 5=1 dipole strength is collected in the giant 
spin-dipole resonance (IVSDR or simply SDR). Several 
investigations on both the experimental and the theoret- 
ical side (see [32| and references therein) have been ded- 
icated to the problem of the GDR-SDR splitting, which 
markedly reflects the dependence of the effective interac- 
tion on the incident energy. A schematic model predicts 
the SDR to be lower in energy than the GDR, on the ba- 
sis of a residual interaction being slightly less repulsive in 
the spin-isospin (5=1, T=l) channel than in the isospin 
(5=0, T=l) one [33j. From general arguments, in a limit 
case one can find a highly collective 1~ state exhausting 
the whole sum rule for the dipole (or for the spin-dipole) 
operator, so that zero transition probability associated 
to the other one is left [Hj|. In practice, this is less likely 
to occur, the two resonances overlap and 1~ states dis- 
playing non vanishing transition strength associated to 
either of the operators can be found to be degenerate or 
lying close in energy (see examples in [Hj|). 

As giant resonances are above threshold, particle emis- 
sion is expected and it can be modelled in TDHF [3(| . A 
discussion about the lifetime of the TDHF wave- functions 
can be found, e.g., in [33], where a comparison between 
the performance of TDHF and the continuum-RPA in 
the matter of escaping widths and the possibility of com- 
paring to experimental information was discussed. 
Ref. [38| presented the effect of two-body correlations on 
the mass dispersion from giant resonances, incorporated 
in the standard TDHF as fluctuation of the one-body ob- 
servable of interest on the basis of the Balian-Veneroni 
approach [39| . which turned out to be important. As a 
matter of fact, although TDHF takes into account the 
coupling to the continuum and those anharmonicities 
that can be captured through one-body operators, exten- 
sions are required to fully describe the spreading width 
of the resonances. Effort has been devoted to formu- 
late extensions of TDHF which go beyond the mean-field 
approximation, capable of modeling coherent collisional 
terms as well as incoherent effects, like those leading to 
2p-2h admixtures into the strength function from non- 
diagonal couplings to low-lying phonons or non collective 



lp-lh bubbles [jof . We recall that couplings of lp-lh to 
high-lying 2p-2h states in spin-isospin channels are ex- 
pected to be favoured by the tensor, a mechanism which, 
spreading strength towards high energies, has been in- 
voked to fully explain the quenching of the Gamow- Teller 
(GTR) and the SDR (see (33[), collective states of par- 
ticular interest for astrophysical processes and particle 
physics. 

If the standard TDHF approach and the extensions re- 
called just above describe a single path in the time do- 
main, stochastic formulations of the theory are available, 
which were shown to be closelyrelated to the Boltzmann- 
Langevin equation (cf. 41 j42l | and references therein. 
Cf. also the discussion in |43j|). 



1. Comments on TDHF and RPA models 

This Section is not aimed at providing a complete re- 
view of similarities and differences between the two ap- 
proaches, their various formulations or the physics in- 
volved. Instead, some (standard) terminology is recalled, 
together with a few elements useful to better follow the 
discussion. Besides the fact that TDHF is a more general 
theory than the RPA, it is useful to remark some techni- 
cal differences when the former is employed in the linear 
limit. 

The RPA framework, firstly formulated for plasma 
physics in 1953 l], involves, besides the construction 
of the static mean-field, the computation of the two- 
body matrix elements of the residual interaction in the 
particle-hole (p-h) channel. In the most microscopic 
approaches, which make little use of ad-hoc adjustable 
parameters, they are respectively defined as the first 
and second derivative of the energy density functional 
(given by ansatz or built on an effective interaction like 
Skyrme) with respect to the particle density fluctua- 
tion [45J. The derivation must be performed from the 
most general expression of the energy density, that is 
built without the symmetry restrictions allowed in spe- 
cial cases (cf. [46j . [25j ; see also 0]). As a matter of fact, 
not only ground-state spin saturation properties (which 
here means ip* (x)<Tiipi(x) = 0) and spatial symme- 
tries are usually broken when exciting the system, but, 
in some cases, the excitation is even mostly driven by 
terms that vanish in the Hartree-Fock approximation. 
Self-consistency, as intended in the sense of the RPA (we 
will not discuss the case of the various types of separable 
RPAs; information can be found in [13]) means that EDF 
contributions are not dropped when proceeding from the 
Hartree-Fock mean-field to the computation of excited 
states. Under this condition, provided the model space 
is complete, it is commonly known that symmetries spon- 
taneously broken by the Hartree-Fock approximation are 
automatically restored and spurious modes, within the 
limits of numerical accuracies, become orthogonal to the 
physical spectrum and degenerate with the ground state 
(zero- modes) [48j . 
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Concerning the TDHF model a full discussion about the 
symmetries breaking should be provided, which is beyond 
the scope of this work. We only add that, as for the RPA 
approaches, numerical inaccuracies can produce spurious 
mixings and alter the strength distribution. However, 
a lack of "consistency" can be produced also in other 
ways, for example from an incomplete implementation of 
the functional at the mean-field level with respect to the 
procedure adopted when fitting its parameters, as well as 
any use which produces lack or overcounting of effects. 

Regarding an effective two-body force with no density- 
dependent coupling constants, as for the tensor terms by 
Skyrme, one can avoid working out the energy density in 
order to derive the mean-field and the two-body residual 
interaction, because, in absence of rearrangement, the 
latter would be identical to the starting force. However, 
one would miss, in this case, the link of the interaction 
terms back to the corresponding energy functional con- 
tribution, to which one usually refers to when comparing 
different calculations. Although it is recognized that the 
globally-fitted values of the Skyrme parameters should 
be interpreted by looking at the functional as a whole, it 
is still interesting to study the role of the various terms 
by setting some coupling constants to zero, or by arti- 
ficially modifying the functional in other ways. If the 
consistency breakings produced in such way is not dra- 
matic, this operation can reveal general features of the 
effective interaction and point to drawbacks and credits 
of one specific parametrization. With such aim, the use 
of Skyrme-like EDFs, that is non two-body interaction 
based EDFs, more or less rich in density-dependent cou- 
pling constants and derivatives, has become quite widely 
employed. 

In the standard TDHF theory, only the mean-field op- 
erator is required in both the static (ground-state) and 
dynamic (residual effects) calculations. The implementa- 
tion of the functional is, in this sense, much less involved 
with respect to the RPA approach and the fullfilment of 
self-consistency less demanding, expecially in comparison 
to RPAs formulated in the coordinate space. The effects 
corresponding to the various terms, coupled one another 
through the densities, are not expected to be simply ad- 
ditive. It occurs in analogy to RPA whenever a different 
matrix is diagonalized once two-body contributions are 
dropped. 

Concerning the microscopic understanding of the re- 
sults, in the widely employed RPA models formulated in 
the configuration space the equation of motion is reduced 
to the diagonalization of a matrix built on a p-h basis; 
the solution immediately provides the microscopic struc- 
ture of the excited states in terms of the simplest excita- 
tions assumed within the theory and this repays the ef- 
fort of building the two-body matrix elements. The RPA 
transition strength distribution is built up by a linear 
combination of the one-body p-h transition amplitudes 
(OBTAs) for the given operator, with weight depending 
on the expansion coefficients resulting from the diago- 
nalization (which reduce to or 1 in the non-interacting 



limit, corresponding to the "unperturbed" response) and 
interference is properly accounted for. These ingredients 
provide together a useful guideline to microscopically in- 
terpret the RPA response (at least in a weak-coupling 
regime). One can understand, for example, which p- 
h wave-function contributes more strongly (in terms of 
OBTA and relative weight) to an eigenstate displaying a 
relevant amount of transition probability for the opera- 
tor of interest. Moreover, one can attest to what extent 
the degree of mixing among the p-h states, the mech- 
anism underpinning the collectivity, where the relative 
importance between the strength of the residual inter- 
action and the unperturbed energies plays a key role, is 
mainly due to one or the other term of the force. The 
reader can refer to one example in [4^ | , which shows how 
the residual Skyrme spin-orbit is able to admix extra p-h 
configurations with opposite spins in the considered low- 
lying RPA state. 

Although less commonly studied, also considered some 
interpretation issues, microscopic information can be ex- 
tracted from TDHF simulations as well, provided the re- 
quested operations are built it. Examples of microscopic 
analyses of TDHF calculations can be found in (50|-[5l|. 
Further considerations on this topic are deferred to the 
future. 

B. Unrestricted Skyrme energy density functional 

TDHF is a well established implementation of the den- 
sity functional approach in nuclear physics § . The start- 
ing point is the static (effective) Hartree-Fock energy den- 
sity H(x) 

E HF = (<S> (X)\H\$ (X)) as = J U(x)dx, (11) 

from which the mean-field can be derived. In our case, in 
the previous equation H is the effective Skyrme Hamil- 
tonian and &o(X) is the Hartree-Fock Slater determi- 
nant. For the sake of simplicity, the kinetic energy will 
not be written in the following, nor the Coulomb term 
(where the exchange is treated within the usual Slater ap- 
proximation). Due to the effective nature of the nuclear 
Hamiltonian and the associated mean-field, the Hartree- 
Fock (Hartree-like) equations share the same philosophy 
of the electronic Kohn-Sham equations [52j . 
Although the Skyrme interaction is constructed as a con- 
tact force, fact that allows the implementation of the 
functional in terms of one-body densities depending on 
a single point of space, finite range effects are included 
through the derivatives of the wave-functions. In fact, the 
standard Skyrme-tensor HF energy density is expressed 
like a sum of terms bilinear in the simple particle or spin 
densities 

=EE #(*»&(*.")U.,' (12) 

i u 

s(*) = EE^( a;, ' w, )^( a: ' a, )< w Vi w )i-=-" 
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where the spin components of the spinor i are repre- 
sented by \oj) (previously indicated by and er are 
the Pauli matrices, and other kinds of densities defined 
on these objects and their derivatives up to the second 
order, which will be introduced further on (see (53j 
for higher order EDFs). Once the static densities are 
replaced by those built on the solutions of the time- 
dependent one-body equations, the total energy density 
gains a time-dependence {T-L(x) —> W.(x,t)), although its 
functional form remains unchanged with respect to the 
static case. 

When only the central and spin-orbit terms of the 
Skyrme force are retained, not all the combinations al- 
lowed by the symmetries, among those depending on the 
second order derivatives, appear in the corresponding 
functional. The tensor force provides a richer structure, 
introducing in the N-N interaction the dependence on the 
relative spin orientation according to 



where r is the relative position of the two interacting 
particles and V{ 2 (r) is the spatial form factor that also 
accounts for the isospin dependence. It is the only (local) 
force able to explain the deuteron quadrupole moment, 
able to transfer to a two-particle wave-function up to two 
units of orbital angular momentum and active only be- 
tween spin triplet states. As a one-body potential, the 
tensor force is able to mix single-particle states that dif- 
fer by one or two units of orbital angular momentum, 
possibly introducing a parity mixing. 
Within a meson exchange model, the form factor of 
Eq. (|13l) is modelled as a Yukawa potential and its zero 
range limit can be parametrized in the Skyrme form [5[- 
1 



v T (l,2) 



4V7 2 (r) 



(13) 



,(<xi -r)(cr 2 -r) 



<J\ ■ <T 2 



T r l T r 

v T (l,2) = — (a x ■ k')(cr 2 ■ k')5(xi - x 2 ) + 5(xi - x 2 ){a x ■ fe)(cr 2 • fe) - — (cti • er 2 ) k' 2 5{x 1 - x 2 ) + 8{x x - x 2 )k 2 

2 1 J o L 

U 



-U((Ti ■ k')8{xi - x 2 ){a 2 ■ k) - y(ci • cr 2 ) k' ■ S(x 1 - x 2 )k 



(14) 



where, as usual, k — ^(Vi — V 2 ) and k' = — 57 (Vi — 
V 2 ) respectively act on the right and on the left. By us- 
ing basic angular momentum algebra, the previous equa- 
tion can be recast in terms of products of tensors of rank 
two in spin and momentum (see eg. fl7^)- 
Although the tensor force is able to act between rela- 
tive L=2 states, where the centrifugal barrier keeps the 
nucleons apart so that the long range attractive part 
of the N-N interaction (r > 3 fm) is probed, there is 
currently no evidence that an intra-medium finite-range 
tensor force performs better than a contact, velocity- 



dependent, parametrization. Some remarks about the 
quality of the zero range approximation for the tensor 
force were proposed by the authors of Ref . 15411 . Recent 
comparisons between Skyrme (SLy5+t of [551 ]) and, in 
particular, the GT2 Gogny force of [56| . can be found 
in ■ Not many finite range effective forces comple- 
mented with the tensor are available and further studies 
are envisaged for the future. 

In the proton-neutron formalism, the contribution to 
the energy density associated to the tensor force (| 14[) . 
once the exchange is taken into account, reads 
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H tens {x) = 2B VS V ■ S n (x)V ■ S p (x) + A vs ]T (V ■ S q (x)f + 2B Jo J°(x)J p ) (x) + A Jo ^(J°(:z)) 2 

i q 

2B. h J n (x) ■ J p {x) + A, h J2(J q ( x )) 2 + 2Bj 2 J n {x)J p (x) + A, h ^(J^*)) 2 + 

g i 
B F [S n (x) ■ F p (x) + S p {x) ■ F n (x)] +A F ^2 S q (x) ■ F q (x) + 

B G [S n (x) ■ G p (x) + S p (x) ■ G n (x)} +A G Y^ S q (x) ■ G q {x) + 

g 

B T [S n (x) ■ T p (x) + S p (x) ■ T n (x)] +A T Y, S q (x) ■ T q (x) + 

i 

B AS {S n (x) ■ AS p (x) + S p (x) ■ AS n (x)) + A AS ^ S q (x) ■ AS q (x). 

I 



(15) 



For a generic term a of the functional, the A a and B a 
coupling constants enter the expression of the proton or 
neutron mean-field, respectively felt by a particle with 
the same and opposite isospin. In terms of the commonly 
employed isoscalar-isovector coefficients C^, they satisfy 
the relations A a = C$ + C? and B a = C^-Cf . The C£ 
coupling constants are related to the Skyrme parameters 
as tabulated in [58], according to the convention for the 
EDF explained below. 

The most evident difference with respect to the standard 
Skyrme functional, based only on the central and spin- 
orbit terms of the force, is the presence of two pseudo- 
vector densities built on the tensor product of the nabla 
operators 

F(x) = ®V) + (V ®V>)]S(x,x')} x=x ,(m) 

G(x) = \{[(V'®V , ) + (y®V)]S{x,x')} m=m ,(17) 



where F M = ±£„[(V 



(V <g> V)] fJ . v S v and 



denotes the usual tensor product between vectors (A Cg) 
B)^ = A^By. The S ■ F and S ■ G terms are connected 
through the (V • S) 2 term according to the equation 



(F(x) + G(x)) ■ S(x) 



-(V-S(x)f 



(18) 



which applies when integration is assumed. The previous 
relation can be used to successfully switch from (TTSJ) to 
the formulation of the EDF of (58|, where the G density 
was not defined and only the S ■ F terms appear, with 
proper weights. In this work we will also work with the 
formulation based on the choice Cy? = 0. More details 



can be found in App.lA] where the derivation of Eq. (fT5|) 
is provided. Indeed, since the terms of the functional are 
not independent from each other, for example also the 
relation [53 



(V-S(x)f 



-S(x)-AS(x)-{VxS(x)) 2 (19) 



holds, many equivalent formulations for the energy den- 
sity can be worked out. 

The S ■ F and S ■ G terms contribute to mantain 
the Galilean invariant properties of the standard Skyrme 
functional when the tensor is made active. In particular, 
they are related to the pseudo-scalar, vector and pseudo- 
tensor densities Jq,J,J_ and the pseudo- vector spin ki- 
netic density T. The latter is defined by 



T{x) = [V- V'S(x,x')} 



x—x' ' 



(20) 



where a generic component depends on V • V'5 M . The 
first three objects are those that represent the trace, the 
antisymmetric and the symmetric parts Q of the pseudo- 
tensor spin current 



J (*) = 2 - [(V - V) ® S(x,x% 



(21) 



which decomposes into 

Jfiv = -^SfivJo + ~ ei ^Ji + J-fivl (22) 



being the Levi-Civita tensor; their square satisfy the 
relations 



/ 2 



J = 



J = 



/IV 

^ ^ J {lis ( J [iv 'Jv{i) 



and simple algebra leads to 
^2 i 



J 2 ; 



(23) 
(24) 

(25) 
(26) 



the usual definition for the scalar product between ten- 
sors AB= Yluv A^B^ has been used. As long as tri- 
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axiality is absent, J = Tr J is zero (in presence of sub- 
scripts, J will be denoted by J^). In spherical sym- 
metric nuclei, the sum over /z, v in Eq. (|25p vanish as 

well [2] , so that one recovers the reduction of J to \ J 2 . 
The Jo , J and J densities enter the central sector of the 
standard Skyrme functional with weights that allow the 
replacement 



Kjrj-i yJ'j 1 ) t~ {J r p *J t ~r ^"v lisp — ^yj* J 'J 1 \£ I J 



for T=0,1, being C^ c = 3C£°' C = 2G 



J\ ,c 

T 



= a 



J°2 , c 



C^' , where c=central, t=tensor). Un- 



( riot A~f Q 

less the tensor is switched on, the two implementations, 

in terms of J |25( or the formulation where Jo, J and 
J are kept separated [58| . are interchangeable. In this 
case, the Galilean invariance only relates the spin current 
terms to those depending on the spin kinetic density T 
(C^' c = —C^?' c ). When the tensor is introduced, the 
following relations among the parameters must hold for 
the Galilean invariance to be respected 



3C£° 
2C^ 2 



7 ■ "T ^{Crp Cji 



a 



c?) 



= 
= 
= 



(28) 



and 3/lOGp 



J ,t 



-2/5C£ ut 



a 



j 2 a 



This is realised for 



the standard Skyrme-tensor functional; in order to prove 
the invariance, it is useful to know that under the local 
Gauge (Galilean) operation Tq : tpi(x) — > e %kx i\)i{x), the 
F and G densities transform according to 

T G {F„(x)) = F^x) + k^k ■ S(x) + k„J (x) + 

+ ^2k v J^{x) (29) 

T G {G^x)) = G^(x)-k^k- S(x)-k tl J (x) + 

-J2 k » J n»( x )- ( 3 °) 

V 

In total, the Galilean invariance imposes constraints on 
eight (six without the tensor) of the fourteen time-odd 
coupling constants of the whole Skyrme-tensor functional 
(in the isospin formalism and in the formulation where 
=0) which, in such a way, remain determined from 
the corresponding even partners. We refer to [25| (notice 
the misprint in Eq. (4.6), where fe 2 has to be replaced 
by k 2 p(r)) and [58] for the left relations concerning the 
Galilean invariance properties and we proceed with some 
other comments. 

The inclusion of the spin current J terms, which, 
as know, impact on the spin-orbit splitting in absence 
of spin saturation (here in the sense that each pair of 
spin-orbit partners is not fully filled up), is required in 
forces that took them into account during the fit. Ne- 
glecting them leads to a breaking of consistency in terms 
of an incomplete implementation of the functional. For 



a more complete discussion refer to 61], where a com- 
parison between predictions from the SLy4 force, fitted 
without that contribution ( "type I" force) , and the SLy5 
set, which included it ("type II"), was considered; ambi- 
guities related to the former type of force were pointed 
out. In particular, it was discussed that the contribution 

from the J terms in computing the excited states can 
be more important than in the ground-state, with the re- 
sult that standard fitting procedures, unable to capture 
some physics through the ground-state observables, can 
lead to poor results when computing collective excita- 
tions. For this reason, suppressing those terms from the 
residual interaction of "type I" forces can lead to an error 
larger than a RPA self-consistency breaking approach, al- 
though (unwanted) case by case considerations with such 
paramctrizations might be required. 
It similarly holds for the time-odd sector of the func- 
tional. The central S ■ T terms have been often sup- 
pressed in the past together with their Galilean partners 

J (cf. e.g. |60(). However, the error associated to the 
suppression of the S ■ T terms in finite systems can be im- 
portant. In [6l| it was shown that the lack of these time- 
odd terms can significantly alter the energy location of 
the GTR, which dominates the landscape in the (charge- 
exchange) 1 + channel. It also turned out that the GTR 
collectivity, which, as known, experimentally absorbs the 
%60 of the Ikeda-Fuji-Fujita sum rule, can be badly af- 
fected by relatively strong (negative) values of the 
coupling constant, a fact which explained the peculiar 
behaviour of the SLy5 parametrization with respect to 
other Skyrme forces. Also, relatively too weak val- 
ues (which depend on the ti and t2 Skyrme parameters) 
would allow a too strong mixing from the S 2 terms (de- 
pending on to and t3, mainly fixed on bulk properties), 
unless the proper balance, or alternative formulations, 
are found. On such basis, one conclusion of [6l| was that 
the spin, velocity-dependent terms cannot simply be ne- 
glected as a rule and the parametrized EDFs rather need 
improved fitting procedures to model the dependence on 
the spin density, possibly including (further) constraints 
from collective properties. 

Among the last developments in the central spin sector, 
extensions of the to-t3 spin-dependent terms were pro- 
posed in [62[ ) i n connection to the spin phase transitions 
that are known to be predicted by effective approaches 
beyond the saturation density (cf. [631 ] and references 
therein). 

The last remark of this Section concerns superfluidity. 
An approximation still quite in use in order to include 
pairing correlations in TDHF consists in performing a 
BCS calculation at the mean-field level and evolving the 
dynamics on top of it. The sums in Eqs. (1121) would con- 
sequently span the larger set of states that defines the 
pairing window, the expression of the ( "normal" ) densi- 
ties gain non trivial occupation factors and the energy 
functional is complemented by terms depending on the 
densities' "anomalous" counterparts 58 1. 
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TABLE I: Central and tensor coupling constant for the J , 

J i J. > S T, S ■ F and J terms, according to the isospin 
formulation for the energy density functional. All the values 
are in MeV-fm 5 . 



TABLE II: At (first line) and Bt (second line) coupling 
constants for the SLy5 force and, within parenthesis, T44 
when the tensor is switched off ("central") or retained. All 
the values are in MeV4"m J . 



SLy5 

central (c) tensor (t) 



T44 

central (c) tensor (t) 



or 


-15.67 


14. 


-59.01 


-24.40 




-64.53 


54. 


-52.03 


20.81 




15.67 


7. 


59.01 


-12.20 




64.53 


27. 


52.03 


10.41 


cf 





-42. 





73.19 


cf 





-162. 





-62.43 


r<J0 


5.22 


23.33 


19.67 


-40.66 


C(° 


21.51 


90. 


17.34 


34.69 


^0 


7.83 


-17.50 


29.50 


30.50 


cl 1 


32.27 


-67.50 


26.01 


-26.01 



A detailed description of nuclear response would require 
a TDHFB approach also to guarantee a proper treatment 
of symmetries. A full 3D-TDHFB is a demanding com- 
putational task. Some of the most advanced works on 
this side are represented by [64]]- 65 1 and, for the Gogny 
case, (66|. A full TDHFB implementation of the Skyrme 



functional, including also the tensor, is not yet available. 
The approximated and not fully self-consistent treatment 
of superfluidity does not affect our conclusions on the p-h 
channel, although an extension which unifies the best ef- 
forts in both the p-h and p-p (particle-particle) channels 
is envisaged for the future. 



1. Parameters, numerical tests and the centre of mass 
correction 

In this work we use the SLy5 force [13] with tensor 
parameters that were perturbatively defined in (55| and 
employed in subsequent applications starting with (68j 
(cf. also [69[ and references therein). As known, the Lyon 
forces, which, in particular, received constraints also from 
microscopic calculations of the neutron matter equation 
of state, were tailored to improve the isovector channel 
of the effective interaction and thus the description of ex- 
otic systems. Like other works, we also employ the force 
T44 of 0], characterized by some different features. For 
example, the central+tensor isovector coupling constant 
C 1 1,c+t of T44 vanishes, so that in the spherical limit one 
is left with the isoscalar tensor contribution only. More- 
over, for all the terms mentioned in the previous Section 
(JC°), J,J,S-T,S-F), the like (C T + Cf) and unlike 
(Cq — Cj) particle tensor contributions to the mean- field 
operator have the same sign, at variance with the SLy5 
case where they are systematically opposite each other 
(see Tab. U. 

The SLy5 and T44 forces have strictly negative central 
C<£ coefficients, a quite common feature of the standard 
Skyrme sets (cf Tab. II in [7(|). They are relatively strong 



central (c) 


tensor (t) 


c+t 


C'o+C'i -80.20 (-111.03) 
Cj-Qf 48.87 (-6.98) 


68.0 (-3.59) 
-40.0 (-45.21) 


-12.20 (-114.62) 
8.87 (-52.19 ) 



for both cases (with a different /Cf ratio), so that the 
effects from some of the newly implemented terms are ex- 
pected to be emphasized with respect to other choices. 
With the inclusion of the tensor, the C^ coupling con- 
stants remain negative; the SLy5 ones receive a large 
reduction (90% and 84% for Cq and Cf respectively), 
while the T44 Cq + Cf and Cq — Cj combinations are 
both strengthened in absolute value (see Tab. [TT| of this 
work) . 

In the pairing sector, an (isovector) volume pairing 
force was implemented. For 120 Sn we set its (neutron- 
neutron) strength to the value Vo(n)=243 MeV-fm 3 ; in 

such a way the SLy5 calculations including the J terms 
well reproduce (<1% of discrepancy) the experimental 
neutron pairing gap (1.32 MeV) obtained with the stan- 
dard mass formula [7l|]. The associated energy cut-off is 
such that, as rule of thumb, only one extra major shell is 
included in the BCS space. In the following, we will re- 
fer to the choice Vo(p)=Vo(n)=243 MeV-fm 3 as "pairing 
set 1" in order to distinguish it from other sets discussed 
later in connection to heavier nuclei. In this work, we do 
not stick to the isospin invariance in the pairing channel 
adopted by some of us in the past [72| and we separately 
fix the proton and neutron pairing strengths. 

The model space in the static calculations is a cubic 
box of 24 3 fm 3 with the commonly employed 1 fm grid 
step. A one-half reduction of the grid spacing would re- 
produce the HF single-particle spectrum with a better 
precision, causing variations from tens to 100-200 keV 
on single-particle energies in 16 O (which do not affect 
the features of the giant resonances we are interested in) , 
but it would importantly enlarge the computational time 
of the dynamical simulation. We consider the standard 
choice of 1 fm to be a good trade-off for our purposes. 
In particular, the convergence in the static was checked, 
as usual, by looking at the stability of the single-particle 
energies. 

The dynamic calculations were performed with reflect- 
ing boundaries over the same grid of the static case, but 
extending the box side parallel to the boost direction 
to 64 fm for the lightest systems. In the other cases, a 
smaller squared box was used to save computational time 
which, together with the other model space parameters, 
had been employed in previous calculations [73| . The box 
size, the simulation time length and the boost strength, 
chosen small enough to stay in the linear regime, min- 
imizing the particle emission but without being dom- 
inated by numerical noise, must be considered mutu- 



10 



ally dependent parameters, properly chosen to avoid un- 
wanted effects from particles reflection. The choice for 
the input parameters can also be tested by checking the 
absence of unphysical fluctuations in the particle num- 
ber or against the stability of the response computed in 
a reduced volume when rescaling the box size. The tests 
were successfully performed by looking at the heaviest 
nuclei. 

The numerical accuracy was also checked by observing 
the free motion of a single ground-state nucleus on the 
grid, the typical test for the Galilean invariance. The 
nucleus in the starting rest condition is boosted by an 
external force and simply put into translation without, 
ideally, internal excitation. 

First of all, one notices that in the simple spin-saturated 
16 O, the contribution to the total energy from the spin 

current J is expected to be zero, the same occurring 
for the energy contributions associated to the time-odd 
terms. Simulating the translation, the former quantity 
turns out to be of the order of 10~ 4 MeV and it oscil- 
lates in time at the 10 -5 MeV level, while the HF single- 
particle energies change by a few tens of keV with respect 

to the "basic" calculation that neglects the J terms. 
The most important error source is due to the choice of 
the grid spacing size: the accuracy improves by one order 
of magnitude after a one-half reduction of this parame- 
ter. In order for the functional to be analytically Galilean 
invariant, the inclusion of the contributions associated to 

the pseudo-tensor J requires, as said, the S ■ T ones too 
and, in presence of the tensor, the S ■ F (S ■ G) terms as 
well. In this system all the time-odd terms vanish with 

accuracy higher than the J terms, leaving no trace on 
the single particle spectrum even at the eV level. 

The last remark concerns the centre of mass (cm.) 
correction to the kinetic Hamiltonian term. In all our 
calculations, the total cm. of the system at rest is lo- 
cated at zero - and it is expected to remain fixed if a 
purely internal excitation is produced. The simple one- 
body approximation was included when fitting the SLy5 
force, so one should keep the correction active in the 
static HF. Concerning the dynamics, it is known [74| - 
[75| that conceptual problems about the definition of the 
mass of the system arise in TDHF, in relation to pro- 
cesses involving more than one fragment (fusion, fission 
and collisions in general). When simulating giant reso- 
nances in TDHF, one needs to deal with the mass disper- 
sion which accompanies the de-excitation process. Some 
other authors computing the GDR [2(| choose to apply 
the centre of mass correction to the mass distribution in 
the whole model space, considering the particle number 
as a constant in the time-space. Other practical recipes 
can be tried; we ran a few test cases in order to com- 
pare a calculation (YN) including the cm. in the static, 
but not in the dynamic, so that some self-consistency 
breaking is introduced, a calculation (NN) that neglects 
it also in the static and the one (YY) accounting for 
the cm. at both stages, with a constant mass number 
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FIG. 1: Isovector dipole (IVD) response in 16 in the time 
domain (upper panels) and after the Fourier transform (lower 
panel) . The comparison among calculations including the cen- 
tre of mass correction everywhere (YY), nowhere (NN) and 
in the mean-field, but not in the dynamical evolution (YN) is 
shown. The vertical scale is arbitrary. For each calculation, 
the expectation value of the lowest order isoscalar dipole op- 
erator is reported as well (dashed line). 



equal to the corrected static one. All these cases con- 
tain approximations, but we wanted to outline which is 
the least worst choice. The outcome from this compar- 
ison is represented in Fig. [TJ where the isovector dipole 
response in 16 O is represented in the time (upper plots) 
and energy domain (lower plot) respectively (the vertical 
scales are arbitrary). The total centre of mass motion 
(obtained from the expectation value of the lowest order 
isoscalar dipole operator Cq-ocb defined analogously to 
O1-01, but with t z replaced by I T , the identity operator 
in the isospin space) is plotted as well (dashed line). In 
the (YN) calculation the presence of a drift of the system 
is recorded and it is possible to notice the spurious con- 
centration of strength at low energies (lower panel). The 
difference between the (NN) and the (YY) calculations 
turns out to be small with the employed parameters in 
the linear regime. 
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III. RESULTS 

As previously recalled, the isovector giant dipole reso- 
nance, which characterizes the experimental cross-section 
more markedly than other vibrational states, is the first 
nuclear collective excitation which was discovered, in 
photonuclear reactions in the 1940's (cf [76]), imagined 
by A. Migdal in 1944 jjjj m connection to the nuclear 
matter polarizability. The microscopic interpretation in 
terms of the (overall repulsive) residual effective inter- 
action, able to explain the collective motion as super- 
position of nucleons undergoing particle-hole excitations, 
was achieved only 20 years later. Since that time, many 
experimental data have become available from a wide 
range of events, together with an extensive theoretical 
investigation. Forty years later than the first observa- 
tion, evidences of finite temperature dipole modes taking 
place in a compound system from a heavy-ion collision, 
were found as well [78} . Some features of the GDR vary 
with temperature (like the size of the deformation and 
the width), while the energy centroids tend to remain 
stable. 

The study of the zero temperature modes, the equiva- 
lent of the Landau zero sound, represent a useful way to 
explore the nuclear structure and also in order to iden- 
tify features of the intra-medium N-N interaction, testing 
how the current models perform. The main purpose of 
the current work is to analyse the effect of the Skyrme- 
tensor functional on (harmonic) dipole states set on top 
of the ground-state, although, in absence of the study of 
the potential energy dependence on the deformation, the 
static mean-field is not guarateed to correspond to the 
global minimum. In particular, in some cases triaxiality 
is forced for explorative purposes. We leave for the fu- 
ture the investigation of finite temperature events, which 
are important in connection to heavy nuclei reactions [79[ 
and in astrophysical environments. 

In this Section, the isovector dipole response from 
TDHF is presented for some benchmark nuclei with a 
mass number ranging from 16 to 238 having, or lead to 
assume, different shapes. As is well known, the IVD re- 
sponse in axially deformed nuclei is characterised by a 
splitting related to the anisotropy of the single particle 
motion: the experimental cross section 

< g x - > (E) = ^<7|i- lO >0E0 + \*\i-,i)(E) (31) 

is a statistical average of the \J n ,K) = |1~, 0), |1~, 1) 
modes, respectively associated to oscillations of protons 
against neutrons along (A"=0) and perpendicularly to 
(K=l) the symmetry axis (see Fig. [2]). Analogously, the 
cross section of giant resonances induced in triaxial nuclei 
is expected to break up into three contributions associ- 
ated to the three characteristic lenghts of the system (k). 
Each distribution can be fitted with a Cauchy-Lorentz 
probability density function 





(Ha (O b 





FIG. 2: Schematic representation of a prolate (a) and an 
oblate (b) nucleus and the relative position of the GDR K=0 
and K=l energy centroids, with a possible outlook due to 
the spreading (scales arbitrary; the height of the peak centred 
around u)b needs to be doubled when comparing to the exper- 
iment). The half-axis are labelled consistently with Eq. (|33l) . 
The symmetry axis (S) points toward the top. 



where T = ET, being T the scale parameter, Eq the peak 
location and ag the corresponding height. From the data, 
one can extract information regarding the type of defor- 
mation (triaxial or prolate, oblate shape) and its size, 
respectively represented by the so-called Hill- Wheeler pa- 
rameters 7 :— arctan v / 2/?22//320 and j3 = \J (3 2 Q + 2/3| 2 , 
built up the dimensionless quadrupole deformation coef- 
ficients /?2 M = if 1<^> ■ The splitting between the vari- 
ous k modes can be described like the shift between the 
energy centroids Ek = ml^/raOfe, computed as the ratio 
between the energy and non energy weighted sum rules 
(EWSR and NEWSR) for each k component (we will 
denote the centroid energy of the whole strength distri- 
bution by E). The splitting is experimentally visible if 
the nucleus is far enough from the spherical shape (,3=0) 
and if it is not masked by other effects [8(| • 

A macroscopic description of the deformation splitting 
is provided by the hydrodynamical model. Depicting the 
nucleus as a spheroid with eccentricity a 2 — b 2 , where the 
semi- axis a is aligned with the symmetry axis, one has, 
in particular, the empirical relation for the quantity = L 

. . Eq 



(see [76] and references therein) 



— = 0.911 v + 0.089, 

u~ 



(33) 



where and u a are the oscillation frequencies of the 
liquid drop induced by an external perturbation along 
the given directions. Table IIIII compares, for 24 Mg, 
28 Si and 238 U, the (SLy5) TDHF energy centroid ratio 



a(E) = a 



(E - E ) 2 + V 2 



(32) 



with the right hand side of Eq. (|33[) based on the HF 
geometry, with the result that the expectation from 
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TABLE III: Comparison between the empirical ratios of the 
dipole K=\ and K=0 energies, based on Eq. (|33p and the nu- 
clear size from the static Hartree-Fock approximation (second 
column), and the TDHF predictions (third column) from the 
SLy5 force (centroids evaluated in the whole energy range). 









A x 


(El/ E())emp 


(Ei/E ) th 


-Mg 


1.385 


1.332 


28 Si 


0.757 


0.758 


238tj 


1.244 


1.240 




10 1 1 ' 1 1 1 ' 1 1 1 

50 100 150 200 250 

A 

FIG. 3: Comparison between predicted (full circles) and em- 
pirical (white circles) GDR energies through different regions 
of mass. The theoretical numbers correspond to the total cen- 
troids E, while the empirical points are from the BF or SJ 
model, as indicated in the main text (the complete empirical 
trends are shown in the inset). 



this empirical evidence is quite well fulfilled. Similarly, 
the energy centroids can be empirically described in 
presence of triaxiality . 

Fig. |3] shows the comparison between the theoretical 
centroids (for all the considered nuclei, despite the de- 
formation) and the predictions from empirical fits. The 
behaviour of the GDR energy when varying the mass 
number is parametrized like hw ~ 31.2^4~ 1 / 3 + 20.6 J 4~ 1/ ' 6 
according to the Berman-Fultz (BF) model, which is 
expected to be more appropriate for spherical light- 
medium nuclei (typically A<100), where the correction 
for the surface (x A^ 1 ^) is more important, or by 
hw ~ 79A- 1 / 3 MeV from the Steinwedel- Jensen (SJ) 
formula (see references in [Hj| , (82[ ) . In particular, for 
each value of the mass number the main plot displays the 
empirical prediction which is closest to the performance 
of the SLy5 force (the full trend from both the models 
is reported in the inset). The BF formula turns out to 
be reproduced better than SJ up to 120 Sn. In any case, 
it is clear that the discrepancy with the result from 
Skyrme is quite high in the light nuclei; although other 
sets can better perform in such cases, the simultaneous 
reproduction of the GDR across different mass regions 
of the periodic table is a long-standing problem [§4| (see 



also [75| and references therein). 

As previously anticipated, in this paper various trunca- 
tions of the SLy5 and T44 Skyrme functional are consid- 
ered, with the purpose of studying how the central spin- 
dependent and tensor terms manifest themselves on the 
nuclear response. Relative effects from the p-h channel 
are therefore discussed. In any case, particular attention 
is reserved to the definition of the pairing force. It should 
be mentioned that the predictions for the percentage of 
total EWSR exhausted in the energy intervals of interest 
are quite sensitive to the model parameters. Higher pre- 
cision calculations are possible in order to provide more 
robust estimates where needed. 

A particular attention is posed on the central J terms 
(cf Sec lII B|) . which had not yet been included in our 
model. In this work we will label by "basic" the calcula- 
tions obtained without the spin current tensor and with 
none of the time-odd (with the exception of the time-odd 
spin-orbit) nor tensor terms. In all the plots display- 
ing the calculated IVD transition strength distribution, 
which is measured in fm 2 unless an artificial smoothing 
is introduced, the scale on the vertical axis has been nor- 
malized to the highest peak produced, for the considered 

nucleus, in the SLy5 calculation including the central J 
terms (labelled by J 2 or simply J 2 . In the plots where 
the K=0 and K=l components are separately plotted, no 
relative factor 2 is included). Comparisons to other the- 
oretical works or experimental data are provided where 
available. 



A. 16 

The first case considered is 16 0, a historically widely 
used benchmark for the TDHF calculations due to the 
low demand of computational complexity. Despite the 
limitations of the mean-field approach and the known 
problems with reproducing the location of the GDR in 
light systems, it still is a good candidate to highlight 
some features of the Skyrme functional. In fact, the resid- 
ual interaction is expected to affect how the strength is 
distributed in light nuclei more markedly than in heav- 
ier systems, where, by a microscopic point of view, more 
single-particle levels are involved in the process and shell- 
dependent effects tend to average out. 

The calculated isovector dipole in 16 O mainly lies be- 
tween 15 and 25 MeV. The upper left panel of Fig. 0] 
shows the response from the SLy5 functional including 

the central J terms (dashed line), in comparison to the 
"basic" calculation defined above (dotted line). They 
produce a shift to higher energy of both the two main 
bumps which characterise the response, peaked at 19.1 
and 22.1 MeV. The most prominent effect is obtained for 
the smaller peak at higher energy, the strength of which 
is reduced by a factor 2 in the more complete calculation. 
An evident change on the strength function appears when 
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TABLE IV: Isovector dipole energy centroids in 16 O for var- 
ious types of calculations not including the tensor. "Basic" 
means that none of the terms specified in the subsequent lines 
have been included. The considered intervals are [17,20] MeV 
(I) and (20,24] MeV (II) for SLy5 ([17,22] MeV and (22,24] 
MeV for the last two rows), [18,21] MeV for T44 ([17,23] 
MeV). The percentage of total EWSR exhausted in each en- 
ergy range is reported within parenthesis. 



SLy5 (I) 



SLy5 (II) 



T44 



basic 

J 2 
J 2 +ST 
J 2 +S 2 
J 2 +ST+S 2 



18.94 (52.3%) 
19.05 (57.0%) 
18.71 (43.6%) 
19.66 (75.2%) 
19.56 (74.1%) 



21.86 (28.2%) 
22.08 (23.7%) 
21.73 (34.8%) 
23.23 (5.20%) 
23.28 (5.90%) 



19.40 (52.1%) 
19.45 (54.4%) 
19.34 (45.6%) 
19.72 (72.3%) 
19.62 (71.6%) 



FIG. 4: Isovector dipole (IVD) strength distribution in 16 
obtained for different truncations of the SLy5 functional. The 
labels list the terms made active in the corresponding calcu- 
lation, in addition to the "basic" version. Details in the main 
text. 
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FIG. 5: Same as Fig. H for the T44 force. 



o2 

including the J Galileian invariant partner in the func- 
tional, that is the S ■ T contribution, which increases the 
response fragmentation for both the Skyrme forces (see 
panels (b)). For the SLy5 case, such calculation pro- 
duces two new bumps of almost equal size, peaked at 18.5 
and 21.5 MeV. This is the case which qualitatively looks 
more similar to the experimental cross section, which, in 
low resolution experiments, displays two main structures 
around 22 and 25 MeV [83[. However, analyses of the 
microscopic structure in such lp-lh limit are required to 
add more information. The panels (d) show the calcu- 
lation obtained when including also the S 2 terms and 
confirm their capability to increase the strength of the 
dipole response, overcoming the effect of the S ■ T terms. 
As a matter of fact, as a consequence of the action of 
the spin squared terms, there is basically one single cen- 



troid with the highest peak at 19.7 MeV, showing an en- 
hanced transition probability and raised in energy with 
respect to the J 2 calculation. The outcomes described 
above are similarly obtained with the T44 force (Fig. [SJ . 
The "aggregating" behaviour of the S 2 ingredient was 
already observed, without the S ■ T terms (panels (c)), 
by Ref. [2(j, where the TDHF dipole response in 16 O and 
Be isotopes was analysed with the SIII force, but no sys- 
tematic calculations were carried on. In general, in a p-h 
view, the spin squared terms, including the rearrange- 
ment due to the density dependence, are able to couple 
spin- flip (fl, l^) , which do not directly contribute to 
the IVD transition amplitudes, and non spin-flip (ff , H) 
particle-hole pairs with configurations of the same family, 
and, to some extent, to also mix them up. By suppressing 
the one-body contribution that depends on the spin 
ladder operators, the effect from the spin squared terms 
on the whole distribution is lost in our TDHF calcula- 
tions and one recovers the two-peaked structure of the 
J c run (the rearrangement has no effect). This finding 
is possibly related to the specific shell structure of 16 O. 
Further analyses will be discussed in the future. 

^2 

The results obtained by adding the tensor J terms to 
the central ones (labelled by J^ +t in the panels (b)) turn 
out to be identical, or displaying a negligible difference, 
to the "basic" version which does not include any of the 
new terms. The inclusion of the tensor S ■ T contribu- 
tions does not produce an appreciable effect too. 
The information about the SLy5 energy centroids for 
the different types of calculations are summarized in Ta- 
ble |TVl together with the percentage of total EWSR cor- 
respondingly exhausted. They are evaluated in the in- 
tervals 17-20 MeV (18-21 MeV for the T44 force), la- 
belled by (I), where the main peak from the J 2 calcu- 
lation is located. The information about the transition 
strength that is left between 20 and 24 MeV (II) is given 
as well. For the calculations including the S 2 terms (last 
two lines), the considered interval is 17-22 (17-23) MeV: 
the strength is collected in a single bump absorbing more 
than the 74% of the total EWSR. In all the calculations, 
a small fraction of strength is concentrated around 10 
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FIG. 6: Time-dependence of the dipole response (upper panel, 
full line) and of the spin-dipole component (central panel, full 
line) described in the main text, for SLy5 in 16 0. The total 
centre of mass displacement (lowest order isoscalar dipole op- 
erator) is plotted as a dotted line in all the three panels. In 
the bottom panel, the same quantity corresponding to the 
worst (YN) case of Fig. [T] is reported. All the quantities are 
directly comparable and measured in fm. 



MeV. The total energy centroid amounts to 20.6 MeV 
(20.8 MeV) for the SLy5 (T44) force. 

Studying the spin response is not the aim of this pa- 
per, nevertheless we tracked in time the expectation value 
of the operator = J^i Dgiic(xi, uji)T z (i), where 

gix(x, oj) = {x)\x\a z (ui). This corresponds to 



<Of-n>(*) 



D 
~2 



A(St l) (a;,i) - Si P \x,t))dx,(M) 



where Si (x,t) — Si P \x,t) is the time-dependent spin- 
isovector density p n = (p„ t - p ni ) - (p pt - p pi ). The 
signal takes place at a scale which is several orders of 
magnitude smaller than the dipole reponse. For some 
truncations of the Skyrme EDF, regular oscillations oc- 
cur from an early stage (FigJHl full line in the central 
panel) and appear to grow up in time. The plot is for 
SLy5 and the same has been found for other cases. We 
verified that within the time length considered, and fur- 
ther beyond, such behaviour, which is not enhanced when 
reducing the box size, is small enough to not influence the 
observed dipole strength. In particular, the oscillations 
are not larger than the spurious total centre of mass mo- 
tion, represented, for the same calculation, by a dotted 
line in all the three panels. This, in turn, is two orders of 
magnitude smaller than the worst case (YN) which was 
shown in Fig. [1] (dashed line in the bottom panel) and 
comparable to the NN and YY cases. Further analyses 
are under way. 
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FIG. 7: Density profile in the x-z plane (symmetry axis along 
the x direction ), describing the late-time evolution (from the 
top-left to the bottom-right) of the instability occurring in a 
J 2 +S 2 calculation in 24 Mg, with the force T44. The cutting 
plane crosses the middle of the nucleus; the darker label 
denotes a density of 0.14 particles/fm 3 , intermediate between 
the contours marked by "3" and "7" . 
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FIG. 8: IVD response from the SLy5 force in 24 Mg, for the 
types of calculations specified by the labels. Only central 
terms are included. 



B. 



The lightest axially deformed system we present is 
24 Mg, where the pairing correlations are negligible. The 
nucleus has a prolate shape, hence one main splitting of 
the response arising from the deformation is expected, 
with the K—0 mode, associated to oscillations along the 
symmetry (major) axis, being less energetic than the 
K=l one. 

The dipole calculations performed with the new terms 
present some important pathologies when using the T44 
parametrization and expecially when the boost is applied 
along the major axis: with the new terms, the system dis- 
plays spurious excitations, enters the non-linear regime 
and is led to fission or explosion. This is the case when 
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TABLE V: Isovector dipole energy centroids in 24 Mg for var- 
ious types of SLy5 calculations not including the tensor, eval- 
uated in the intervals [0,18) MeV (I), [18,23) MeV (II) and in 
the larger [18,26] MeV window (III). The exhausted percent- 
age of the total EWSR is listed within parenthesis. 



SLy5 


I 


II 


III 


J 2 


15.81 (30.0%) 


20.14 (31.5%) 


21.78 (56.5%) 


J 2 +S 2 


15.85 (29.9%) 


20.74 (44.2%) 


21.42 (56.8%) 


J 2 +ST+S 2 


15.79 (30.1%) 


20.58 (42.0%) 


21.37 (55.5%) 



adding the tensor to the central J terms: the dipole os- 
cillations are dramatically enhanced and the system splits 
up like as result of an inelastic collision. Rotational ef- 
fects can arise as well, like in the J 2 +S 2 model with no 
tensor (Fig. [7]), where the dynamic is similar to a fusion 
state after a non-central collision. With the J 2 +ST+S 2 
model, octupole deformations appear on top of the dipole 
oscillations. When the boost is applied along the perpen- 
dicular direction, the pathologies due to the inclusion of 
the S 2 terms are largely reduced when the S ■ T terms 
are added as well. 

When using the SLy5 force, it is necessary to run the 
simulation longer (at least for a time doubled with re- 
spect to the J 2 +S 2 case), to face the same kind of in- 
stabilities, which thus manifest outside the time window 
considered in this work, far enough to not affect the re- 
sults. These are shown in Fig. [8j Non zero signal is 
mainly visible between 13 and 30 MeV, with a narrow 
peak associated to the K=0 mode located at 15.8 MeV 
(practically unchanged when switching on/off the tensor 
and time-odd terms of the functional), and a broader 
bump associated to the second characteristic length of 
the system (K=l), with a narrower structure in the 18- 
23 interval (II). The central J 2 model performs similarly 
to the "basic" one. For the runs without tensor, the plot 
presents a pattern similar to Oxygen: the S 2 terms are 
able to reduce the spreading of the dipole response with 
respect to the central J 2 calculation (dashed line), while 
the S ■ T terms work in the opposite direction (full cen- 
tral calculation J 2 +ST+S 2 represented by a full line). In 
the inset, the J+ST calculation is separately drawn. 
All the information about the energy centroids computed 
in the intervals of interest are listed in Table [Vj together 
with the corresponding fraction of exhausted EWSR (see 
the caption for details). The calculation obtained by 

adding only the J terms to the "basic" places in the 
interval II about a 10% of energy-weighted strength less 
than the other cases, because the K—l mode is frag- 
mented at higher energies. The total centroid amounts 
to 20.2 MeV for all the considered truncations of the SLy5 
functional. 

The inclusion of the time-odd terms reduces the over- 
all deformation splitting by 300-400 keV, which glob- 
ally amounts to about one half of what is predicted by 
the Gogny-based QRPA of Peru and Goutte [H|. This 
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FIG. 9: K=l (left panels) and K=Q (right panels) compo- 
nents of the IVD strength distribution in 28 Si from the SLy5 
(top) and T44 (bottom) force. The comparison between the 
"basic" and the (central) J 2 calculations is shown. The "J 2 

on basic" run has been obtained by suppressing the J terms 
in the ground-state. The arrow indicates the height of the 
main peak in the latter case. 



change is only due to the residual effects in the p-h chan- 
nel, since the time-odd terms vanish in the ground-state 
and do not alter the single-particle spectrum. Our out- 
come is similar to the SkM* [86[ (Q)RPA result presented 
by |87| (same deformation parameters and r.m.s. are ob- 
tained in the ground-state). 



28 Si is a system comparable to 24 Mg in mass, but with 
different deformation properties. Both the SLy5 and T44 
forces predict an oblate shape (7=60°), with a j3 param- 
eter that, with the second force, is about 30% weaker 

when the central J contribution is not included (/3=0.19 
in the "basic", versus /3=0.26 of the J 2 case); it increases 
by a further ^10% when the tensor is switched on. The 
variations are smaller in the SLy5 case (/?=0.28±0.01 for 

the three calculations), with the J still enhancing the 
deformation. 

Fig. [9] shows how the linear response changes between 
the J 2 run and the "basic" version for the SLy5 (upper 
panels) and T44 force (bottom panels). The K=l and 
K=0 modes are separately computed. It is possible to 
notice that the peak around 22 MeV receives contribu- 
tion from both the K=0 mode and the Landau spreading 
of the K=l component. For the T44 case, the spin cur- 
rent terms enlarges the splitting by 1.4 MeV, since the 
K=l centroid lowers from 20.5 to 19.5 MeV, while the 
K=0 one increases by 400 keV, as reported in Tab. ED 
as a consequence, the ratio between the centroid energies 
decreases by 6% (third column). The empirical estimates 
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TABLE VI: K=l and K=0 energy centroids from T44 in 28 Si, 
evaluated in the whole energy range, in the two self-consistent 

calculations where the (central) J terms are neglected (ba- 
sic) or included in both the static and the dynamic evolution 
(J 2 ) and in the basic+J 2 case explained in the main text. 
The comparison with the empirical ratio of Eq. (|33[) is given 
as well. 



T44 


E 1 [MeV] 


E [MeV] 


(Ei/E )th 


(Ei/ Eo) emp 


basic 


20.45 


24.38 


0.839 


0.837 


J 2 


19.45 


24.79 


0.785 


0.780 


basic+J 2 


19.99 


23.82 


0.839 


0.837 
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FIG. 10: IVD response in 28 Si for various types of SLy5 cal- 
culations. The white inset shows the (rescaled) profile of the 
response perpendicularly to the symmetry axis (K=l). 



TABLE VII: Isovector dipole energy centroids in 28 Si for the 
SLy5 force evaluated below (I) and above (II) 20 MeV, for 
the calculations listed in the first column (no tensor). 



SLy5 


I 


II 


J 2 


17.13 


25.16 


J 2 +ST 


17.02 


24.85 


J 2 +S 2 


17.38 


25.25 


J 2 +ST+S 2 


17.40 


26.02 





1 1 


1 1 1 1 1 

1 <J > , 




28 Si 








— (r+ST+s") c+l 






"j (f+ST+S 2 +SF) c+i 


0.6 












£ 0.4 


— 1 






1 




0.2 


f 






" , J 


^^^V-^— T44 










I L L^J I I T— I i I 

10 20 30 40 50 
E [MeV] 



FIG. 11: IVD respose in 28 Si for the T44 calculation listed in 
the legend. 



of Eq. (I33|) are listed in the last column for comparison. 
The effect from the SLy5 force is less pronounced, but 
similar. 

As said, a useful tool of investigation of how the en- 
ergy functional reproduces the nuclear response consists 
in performing not self-consistent calculations where the 
underlying mean-field is obtained from a functional sim- 
pler than what used in the dynamic evolution. In such 
a way one can disentangle effects arising from the static 
mean-field, including the change produced on the intrin- 
sic deformation. By doing this in the T44 case, it is 
possible to see that the 1.4 MeV variation of the shift 
between the K=0 and K=l modes can be explained in 
terms of the change produced on the ground-state, which, 

as said, is made more oblated by the J . In fact, when 
turning the spin current on in a dynamical calculation 
performed on top of the "basic" mean-field (basic+J 2 .), 
one obtains a profile similar to the fully self-consistent J 2 , 
run (see the dashed lines in Fig. [9] for T44), but the gain 
in the energy splitting is practically lost and one recovers 
exactly the same energy ratio of the pure "basic" case (a 
small fraction of strength can be noticed at low energies, 
below 10 MeV, due to the self-consistency breaking). 
From Fig.[TUl showing the total IVD strength from the 



SLy5 force, one can identify, once again, the opposite be- 
haviour of the central S ■ T (upper-right panel) and S 2 
terms (lower-left panel) in redistributing the transition 
strength. Table I VIII displays the energy centroids com- 
puted below (interval I) and above (II) 20 MeV. With 
the S 2 ones, strength is pushed at higher energies. The 
S ■ T terms fragment the response, with the result that 
the main peak is lowered by 100 keV and transition prob- 
ability is now visible between 18.0 and 21.5 MeV. 
Concerning the tensor, instabilities arise when including 
the S ■ F terms. In the case of the SLy5 force, these ap- 
pear at the middle of the time length we have typically 
assumed. We still attempted to compare with runs where 
the coupling constants are zero by halving the simula- 
tion time. This operation might lead to loose relevant in- 
formation, although it turns out that, when repeating the 
same operation with the cases listed in Tab I VIII the rela- 
tive differences between the calculations still hold. From 
the fourth plot of Fig. [TU1 it seems the S ■ F terms pro- 
duce an overall attractive effect. The same comparison 
between T44 calculations including the tensor is shown in 
Fig. QT] The S ■ F terms shift down the highest peak by 
a further ^300 keV and also move strength (if=0 mode) 
to higher energies. In any case, the role of these tensor 
terms should be further investigated. 
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FIG. 12: The right hand panel shows the dependence of the 
IVD response in 120 Sn on the pairing set "1" and "2" defined 
in the main text, when using the T44 force. The left hand 

++2 

panel shows the effect of the inclusion of the tensor J terms 
(dashed line), in comparison to a calculation including only 
the central counterpart (full line) and when none of them is 
taken into account (dotted line). 



FIG. 13: IVD response in 120 Sn from the SLy5 force (no ten- 
sor) ; the three panels show the change produced when includ- 
es 

ing the J terms (left panel, full line), when also adding the 
Galileain invariant partners S T (central plot, dashed-dotted 
line) and when switching on the S 2 ones as well (right panel, 
dashed line). 



D. 120 Sn 

The next case under study is the spherical 120 Sn, 
where, at variance with the previous cases, pairing cor- 
relations are important and are modelled with the pair- 
ing parameters (set 1) introduced before. Superfluidity 

<->2 

opens the channel to the J tensor terms also in the 
ground-state, preventing the N=70, Z=50 core from be- 
ing spin-saturated, although the system remains time- 
reversal invariant. Consequently, we still expect no con- 
tribution from the time-odd terms in the mean-field. 
While checking this, we noticed that the inclusion of the 
S ■ T terms may help the convergence in the static cal- 
culation when the squared spin density is kept included. 

When the J terms are switched on, the pairing gap 
changes by 1-2%, which is far from producing a rele- 
vant change on the main properties of giant resonances. 
When using the optimum SLy5 pairing set in connection 
to the T44 force, a pairing gap of 1.84 MeV is obtained. 
However, even a change of around +20% of the pairing 
strength ("pairing set 2"), which significantly alters the 
pairing gaps, would produce small variations on the re- 
sponse (see Fig. [T^l left panel). There is consequently no 
need to further tune the pairing either when modifying 
the functional for a given force or when changing from 
SLy5 to T44. 

The SLy5 strength distribution is characterized by two 
main peaks (see Fig. The predictions for the energy 
centroid (-15.4 MeV, calculated in the 13.0-18.5 Mev 
interval) match the experimental value reported in [83[ 
(cross-section data fitted in the 13-18 MeV range). The 
change produced by the new terms in this superfluid 



A=120 nucleus is less evident than in the lighter sys- 
tems. However, one can still notice that the S ■ T terms, 
in the central plot, (slightly) reduce the peaks' height ra- 
tio and produce an attractive effect on the centroids, at 
variance with what obtained when adding the S 2 terms 
(right panel). All the calculations modify the low energy 
tail of the strengh function (below 13 MeV), although 
these outcomes must be considered carefully, since the 
low lying states, besides being sensitive to higher order 
correlations and the pairing, can be particularly affected 
by numerical artifacts from the Fourier transform and 
spuriousities in the wave-functions. 

When the tensor J and S ■ T are made active, no noti- 
cable change in the strength function is introduced when 
SLy5 is used. This is not the case with the T44 force, 
where the J^ +t calculation is somewhat different (see 
Fig. [121 right panel). We finally notice that the two- 
peak structure always visible with the SLy5 force is lost 
when using the T44 set in the simplest "basic" calcula- 
tion, where the distribution is shifted to lower energies, 

but it appears as soon as the J contribution is added. 



E. 190 W 

In triaxial nuclei the deformation leads to a more com- 
plex nuclear response when the system is excited. The 
intrinsic nuclear shape, which impacts on the dynamical 
response, results from the balance of different features of 
the intra-medium N-N force: on one side, various trun- 
cations of the Skyrme functional in the p-h channel can 
affect the deformation properties in the ground-state; on 
the other hand, (monopole) pairing correlations tend to 



18 




steps steps 



FIG. 14: Evolution of the deformation parameters f3 and 7 in 
190 W during the first stage of the iteration process in Hartree- 
Fock, for SLy5. Different lines correspond to different runs 
complemented by the pairing force 1, as explained in the leg- 
end, ft is adimensional, 7 is measured in degrees. 



TABLE VIII: Isoscalar+isovector Hartree-Fock energies (in 
MeV) from the central (J (0) ) 2 , (J) 2 and (J) 2 terms of the T44 
functional, in the prolate (7 = 0°) and triaxial (7 = 39.3°) 
minima found in 190 W. 



T44 


Ai?./ 


AE, h 


A£j 2 


7 = 0° 


x 1(T V 


11.951 


0.157 


7 = 39.3° 


x 1(T 5 


9.671 


0.105 



drive the nuclear system towards a spherical shape and 
alter the deformation parameters as well, so that the 
proper focus on the pairing force must be posed when 
studying intrinsically deformed systems where superflu- 
idity plays a role. 

The measurement of the E 4 + / E 2 + energy ratio, which is 
a good signature of shape transition, suggests that 190 W 
behaves like a triaxial rotor [88[. When using the op- 
timum 120 Sn pairing force (set 1) in connection to the 
SLy5 functional, a triaxial shape is obtained from our 
calculations. The two panels in Fig. H~4l show (an extract 
of) the evolution of the (3 (up) and 7 (down) parameters 
in 190 W with the number of iterations with the SLy5 
functional, for a "basic" calculation and when the spin 
current pseudo-tensor is included (with and without ten- 

■H-2 

sor). This is an example where the inclusion of the J 
terms quicken the convergence. In particular, a 14% de- 
crease of the 7 value (21.8°) with respect to the "basic" 
calculation (25.3°) is produced. A SLy4 calculation ob- 

tained with no J , with the same starting conditions, 
shows a 7 parameter (7=23.2°) close to the J 2 +t SLy5 
one (7=23.7°). The (3 parameter does not change. 
The type of truncation of the energy density functional 
can not only induce a variation of the deformation pa- 
rameters for a given (local or global) minimum, but also 
reverse the relative location in energy of different min- 

<->2 

ima. For example, when the J terms are neglected in 
the functional based on the T44+pairing 1 parameters, 
it is possible to identify one axially deformed and one 
triaxial solution, with an energy difference of 0.91 MeV, 
where the prolate configuration is the most bounded. By 
switching those (central) terms on, the situation is re- 
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FIG. 15: Evolution of the deformation parameters of W 
on the pairing force, for the T44 Skyrme set with the iter- 
ation step number. The pairing sets used in this work are 
considered. 



versed, with the triaxial configuration becoming deeper 
by 1.42 MeV than the prolate shape (with small varia- 
tions of the /3 and 7 parameters, respectively amounting 
to ~1% and <~3%). This extra binding is mainly due 
to the fact that the (positive) energy contribution from 

the J terms in the triaxial case (9.78 MeV) is 2.3 MeV 
weaker than in the axial state (see Tab. IVIlTl) . 

Turning on the time-odd terms, the Kramer's degener- 
acy is expected to be preserved. Their inclusion tends to 
make the convergence more difficult and the time-reversal 
invariance is easily broken; in particular, the use of the 
T44 force required a particular effort to find suitable ini- 
tial conditions. In some cases the problems are empha- 
sised when the S ■ T terms and the S 2 ones are included 
separately. 

We encountered a similar situation when changing the 
pairing parameters, including the values of 288.5 and 
298.8 MeV-fm 3 for the neutron and proton pairing 
strengths respectively ("pairing set 2"). These are opti- 
mal parameters in connection to the SLy5 force, as they 
provide pairing gaps A n p <E 0.7-0.8 MeV, much closer to 
the empirical estimates (0.859 and 0.740 MeV) than the 
quite small pairing gaps obtained with the pairing set 1 
(A, £ 0.2-0.4 MeV, q=n,p). 

We notice interesting variations on the deformation prop- 
erties when changing the pairing force. Fig. [15] shows 
the outcome from "basic" T44 calculations based on the 
pairing sets considered in this work. If the choice of the 
pairing set 2 leads, in this case, to a prolate axially de- 
formed shape, by reinforcing the proton pairing strength 
by -15%, that is up to 344.00 MeV-fm 3 ("pairing set 3"), 
one triaxial minimum with 7=42° is constructed. Revers- 
ing such proton and neutron pairing strengths ("pairing 
set 3a" ) would force the system to converge to an oblate 
shape (dotted line; notice the spike around the 100th 
step). 

Concerning the projection of the linear response in the 
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FIG. 16: IVD transition strength distribution in 190 W on 
top of a triaxial ground-state. The x (main axis), y and z 
components are separately plotted for some truncations of 
the SLy5 Skyrme functional. 
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FIG. 17: IVD strength distribution in 190 W from the T44 
force, in a "basic" calculation (left panel) and when includ- 
ing the tensor at the J 2 (middle) and J 2 +S 2 +ST level (right 
panel) . 



IVD channel, in a triaxial case we expect a three-modal 
response associated to the characteristic lengths of the 
nucleus. Fig. [T7j] shows the transition strength distribu- 
tion in 190 W from the SLy5 force, along the x, y and z 
Cartesian axis and for some truncations of the Skyrme 
functional. Each mode appears as a bi-modal function, 
where the highest peak is respectively located at ~12, 
13 and 14 MeV and the weakest structure is around 16 
MeV. The latter bump becomes increasingly more impor- 
tant from x to z while the main peak reduces in height by 
^40%. Concerning the various functional terms, the S-T 
ones (Fig. [T6l upper panels) produce an overall (small) 
attraction on the y and z components. The inclusion of 
the S 2 terms (Fig. [T5] bottom panels) tends, instead, to 
remove the bi-modal structure in all the components and 
the extra mixing places the (new) centroids at higher en- 
ergy. 

The comparison of the total strength distribution among 
three types of calculations including the tensor with the 
T44 force is shown in Fig. [T7J The effect of the cen- 
tral+tensor J 2 and J 2 +S 2 +ST runs consist of a broad- 
ening of the response with respect to the "basic" case. 
Although the inclusion of the tensor can help the nu- 
merics in some cases, no stable calculation including the 
tensor S ■ F terms is currently available in this nucleus. 

F. 178 Os 

Another system under study is 178 0s, which has been 
considered a good candidate for shape cohexistance (see 
[soj ] and references therein). As for 190 W, we exploit the 
7-softness property to perform theoretical investigations, 
by forcing triaxiality in some cases. 
In this system, both the empirical proton and neu- 
tron pairing gaps belong to the interval 0.91-0.93 MeV. 
Among the various pairing sets considered, the force 2 
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FIG. 18: Same as Fig.HU for 178 0s. 

still appears the most suitable, although the convergence 
is not trivially ensured, expecially for T44. The SLy5 
force behaves better also in this system, in terms of con- 
vergence. The HF calculation easily lead to a triaxial 
minimum; Fig. [18] shows that with the pairing force 1 
(to which pairing gaps A„=0.7 and A p =0.4 MeV corre- 
spond), the 7 parameter changes by ~13% when turning 
the spin current terms on. 

When performing TDHF calculations on top of the 
available mean- fields, we find outcomes similar to the 
triaxial case of 190 W (see Fig. [TO]): the central S ■ T 
contribution provides an attractive effect on the y and z 
components and increases the spreading at variance with 
the S 2 terms. Along the main axis, with the S ■ T terms 
the signal relaxes down into the tails, which become fat- 
ter while the height of the main peak is reduced by ~18%. 
Fig. [20] shows the comparison between a purely central 
J 2 calculation and the one including the tensor in a full 
way, with an evident difference in terms of fragmentation, 
especially on the y component (the total EWSR is still 
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FIG. 19: Same as Fig.HH in 8 0s. 
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FIG. 20: Projection of the SLy5 IVD response along the three 
axis of 178 Os, when the tensor is fully included in comparison 
to a purely central J 2 calculation. The arrow indicates the 
height of the main peak in the former calculation. 






















FIG. 21: The same as Fig. [7] for a SLy5 calculation in 178 0s 
including the central+tensor J 2 +S 2 +ST+SF terms, at a late 
stage of the dynamics. 
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FIG. 22: IVD components in 238 U for SLy5. 



preserved). Running the simulation longer, instabilities 
occur at the end: rotational currents set up and finally 
the system expands itself in the whole model space (see 
Fig. El). 



G. 



The heaviest system considered in this work is the ax- 
ially deformed 238 U (/3=0.27 with the SLy5 force). The 
inclusion of the tensor S ■ F contribution leads to a stable 
result, which is shown in Fig. [22] (full line) in comparison 
to a simpler J 2 calculation (dashed line). The K=0 and 
K=l modes are separately plotted at the right and left 
hand side respectively. 

The resulting smoothed total strength distribution is 
shown in Fig. [521 The centroids of the K=0 and K=l 
modes are located at 11.0 and 13.5 MeV respectively, not 
far from 10.92 and 13.98 MeV reported in [83j]. 
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FIG. 23: IVD response in 238 U from SLy5 with a full inclusion 
of the tensor force. 
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After the smoothing, which simulates the broadening ef- 
fects due to higher order effects not accounted for in 
TDHF, the peaks just above 15 MeV are not resolved any 
longer. Such level of detail is dependent on the employed 
box (discretization) and other possible numerical arti- 
facts. However, one cannot exclude that, when using dif- 
ferent parameters, the full treatment of the tensor might 
more markedly affect the left hand side of the strength 
distribution, which, in any case, appeared quite sensitive 
to the Skyrme set in calculations where the tensor was 
not included. For example, in Ref. [90| . where a simpler 
functional based on the SkI3 set was employed, and in 
the separable RPA calculations of [9l|, based on the pa- 
rameters from SLy6, a shoulder at the left hand side of 
the transition strength function was found. Finally, we 
recall that the modeling of higher order correlations is 
expected to alter the tails of the strength function. Such 
extensions with the full Skyrme-tensor functional are un- 
doubtedly envisaged for the future. 



IV. CONCLUSIONS 

This work presents, for the first time, the full im- 
plementation of the Skyrme energy density functional 
(EDF), complemented with the original formulation of 
the Skyrme-tensor force, in the three-dimensional time- 
dependent Hartree-Fock (3D-TDHF), with no symmetry 
restrictions. To our knowledge, this is the most complete 
development of TDHF concerning the p-h channel, where 
full self-consistency is achieved between Hartree-Fock 
and the dynamics (residual interaction). The derivation 
of the Skyrme-tensor Hartree-Fock energy is provided as 
an extension of [25[ . 

The model has been applied in the linear limit, where 
the random-phase approximation (RPA) is recovered, to 
simulate the nuclear response in the 1~ channel, in- 
duced by the isovector dipole (IVD) operator. Comments 
about the methodology have been outlined and some self- 
consistency issues discussed. 

Various benchmark nuclei with a mass number ranging 
from 16 to 238 have been considered and triaxial calcu- 
lations have been presented. The overall results satisfy 
the empirical trend from the liquid drop model, although, 
as expected, deviations are found in the light region of 
mass, where the giant dipole resonance (GDR) energy is 
underestimated. More importantly, the gross features of 
the strength distribution from TDHF have been found 
in agreement with the experimental data or other theo- 
retical works for all the systems where information are 
available. 

The effects produced by each spin-dependent term of the 
p-h sector of the functional have been separately dis- 
cussed. These are more visible for lighter nuclei where, 
in the RPA picture, fewer p-h configurations are involved 
in the response and Landau fragmentation is more pro- 
nounced. However, systematic outcomes are found, for 
both the considered Skyrme forces (SLy5 and T44), from 



16 to 238 U: the GDR transition strength, which re- 
ceives contributions from one-body transition densities 
(OBTDs) between parallel spin orbitals, is appreciably 
sensitive to the central spin-dependent part of the resid- 
ual interaction. In particular, a behaviour similar to cal- 
culations in other (charge-exchange) channels sensitive to 
the spin (-isospin) residual force has been noticed. This is 
the case of the Gamow- Teller resonance, modelled within 
a RPA formulated in the configuration space 61], the 
transition strength of which receives contribution from 
both spin-flip and non spin-flip OBTDs. The S 2 terms 
tend to reinforce the residual interaction in the isovec- 
tor channel, providing a repulsive effect on average and 
enhancing the strength of the resonance markedly in the 
considered cases. This dominant mechanism is countered 
by the S ■ T terms, which, although less important in 
the IVD channel than what was found in the (charge- 
exchange) 1 + one, produces more fragmentation in the 
GR region. 

On the recalled basis, one would conclude that the 
S ■ T terms, which are spin and velocity-dependent, can 
be useful to balance the S 2 contribution of the current 
fuctionals. The former have been often suppressed in the 
past, together with their Galilean partners, the time-even 

J terms. In Ref. [92J, the removal of the S ■ T terms 
was found to improve the performance of the Skyrme 
EDF in infinite matter; however, as acknowledged by the 
same authors, this suppression would completly remove 
the Gj =0 and Gf =1 Landau parameters of the standard 
Skyrme functional, besides spoiling, with more or less 
manifested effects, the Gq =0 , Gq =1 ones, which depend 
on both the to, ta and ti, t2 parameters, respectively 
through the and coupling constant. 
Besides influencing the spin-orbit splittings, the spin cur- 

rent terms ( J ) are of particular interest in connection 
to tensor studies, as they allow the tensor to become ac- 
tive even under time-reversal invariant conditions. For 
this reason, they were the first tensor terms to be consid- 
ered in the past @, as well as in the more recent history 

0, M- 

<->-2 

In the ground-state, the inclusion of the J terms can 
markedly influence the deformation properties. In par- 
ticular, variations of the deformation parameters up to 
~30-40% have been recorded in 28 Si, where the splitting 
between the K=0, 1 modes accordingly changes by ~1.5 
MeV, and 178 Os. In 190 W, they turned out to drive the 
system to a triaxial shape. 

The inclusion of the tensor at the J and S ■ T level does 
not introduce important changes on the IVD response in 
the considered cases. A small effect is produced by the 
S-F terms, but higher precision calculations are required. 

A particular attention has been paid to the pairing cor- 
relations, which affect the intrinsic deformation in terms 
of both the j3 and 7 parameters. This has to be taken 
into consideration when discussing the predictive power 
of EDF-based calculations. More precisely, the impact 
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of the smearing of the Fermi surface produced by a zero 
range (monopole) pairing force has been discussed: dif- 
ferent values of the pairing parameters have been consid- 
ered, tuned against experimental masses or exploited as 
additional degrees of freedom to force the intrinsic nu- 
clear shape. The treatment of the p-p channel does not 
alter our conclusions about the p-h channel, where rel- 
ative effects are mainly considered, and do not prevent 
the comparison with the experiment. 

The occurrence of instabilities from the Skyrme func- 
tional has attracted particular interest in recent times 
(cf. [93[ and references therein). In our calculations, var- 
ious types of instabilities from the spin-density dependent 
terms of the SLy5 or T44 functional have been found, 
which are enhanced in presence of derivatives, depend 
on the way the EDF terms couple one another and mix 
with the zero modes. They appear similarly to the pos- 
sible cases described in [93| ■ This is especially true when 
the T44 force is employed, while the SLy5 set behaves 
more regularly. A peculiar case is represented by the 
performance of T44 in 24 Mg, where different truncations 
of the functional allow the spuriousities to dominate the 
scene in a relatively short time. This might be related 
to possible softness properties of this nucleus, which are 
predicted by some calculations [87[. The T44 force, more- 
over, turned out to hinder the convergence in the static 
Hartree-Fock quite easily in our calculations, making the 
numerics more difficult. 

We notice that the S 2 and the S ■ T terms can pro- 
duce an unphysical time-reversal breaking in both the 
ground-state and the long term dynamics of the heavy 
7-soft nuclei we considered. The inclusion of both those 
EDF contributions turned out to reduce the appearance 
of unwanted effects in some cases, with respect to the 
situation where only one of them is included. It similarly 
happens when the tensor is added to the central S ■ T 
terms, due to cancellations. The terms depending on the 
laplacian and on the divergence of the spin density have 
not been analysed in this work. No instabilities related to 
the laplacian of the particle density or to the other time- 
even contributions have risen up, but more focused and 
systematic analyses in both sectors can be accomplished. 

Unphysical dissipation effects in TDHF, when time- 
odd terms from the spin-orbit sector of the functional are 
suppressed, were observed in [36] and divergences taking 
place when running the simulation for a long time, in 
relation to couplings with zero-modes, were mentioned 
in [20] . A more detailed study of such phenomena is en- 
visaged also in view of computing the spin modes and 
we will deal with this subject in a separate work. The 
3D-TDHF could represent a particularly sensitive test- 
ing tool for Skyrme(-like) energy density functionals in 
finite systems and complement information from infinite 
matter analyses (94|-[95j]. Whereas the problem of insta- 
bilities might represent a limitation and demands for a 
cure, also to not be obliged to drop terms and possibly 
introduce further spuriousities, useful information about 
the nuclear structure can be accessed through correla- 



tions associated to the symmetries breaking. 

Changing the functional by dropping terms according 
to the one's needs is useful for explorative purposes, the 
ultimate goal of this line of research still being the real- 
ization of a reliable functional, ideally flexible enough to 
adapt to disparate conditions. Finding the proper bal- 
ance between the various spin-dependent terms would 
improve the predictive power of the current formulations. 
As is known, one conclusion of Ref. was that the avail- 
able parametrizations for the tensor force are not able to 
adjust the drawbacks of the central and spin-orbit part of 
the Skyrme functional, although it can offer extra degrees 
of freedom to fix the velocity-dependent terms. Interest 
is still driving the search for suitable constraints sensitive 
to the spin-dependent terms from excited states and/or 
ground-state properties. Clearly, in order to ensure pre- 
dictability on collective excitations in the RPA or equiv- 
alent approaches that strongly depend on the underlying 
shell-structure, coherent efforts are required to improve 
the description of the involved single-particle degrees of 
freedom. 

The GDR energy is known to be correlated to the sym- 
metry energy in dependence on the TRK enhancement 
factor, expressed in terms of the ti and t2 parameters 
through the effective mass (96[, which is already em- 
ployed in fitting procedures (see e.g. [67|). Input from the 
deformation splitting could also be taken into account. 

Concerning further applications, we saw that the spin- 
dependent terms of the Skyrme functional can affect the 
tails of the nuclear response. The fraction of dipole 
strength concentrated in the low energy tail of the GDR, 
around the threshold for particle emission, impacts on 
the photodisintegration rates, of particular interest for 
nucleosynthesis (see e.g. [97j). In such region, also inter- 
ested by pigmy resonances, pairing correlations and the 
coupling to low- lying phonon can play a role [98j. The 
investigation of such energy window by means of TDHF 
simulations, even in the linear approximation, is delicate, 
because proper attention is required to remove artifacts 
related to the Fourier transform and possible unphysical 
effects from the coupling to zero modes. In any case, the 
treatment of anharmonicities with the full Skyrme-tensor 
functional, in particular the extensions to account for the 
coupling to high lying 2p-2h states, which is expected to 
be strengthened by the tensor force, would be an inter- 
esting step to be undertaken. 

Finally, studying particle emission from collective excita- 
tions allows one to access information about the involved 
single-particle degrees of freedom as well as the mech- 
anisms responsible for their cooperation, of interest for 
both the (interdependent) nuclear structure and nuclear 
interaction modeling. 
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Appendix A: The unrestricted Skyrme-tensor Hartree-Fock energy density functional 



In this Appendix, the derivation of the Skyrme energy density functional provided in [25| is extended in order to 
include the tensor contribution, under the general hypothesis of no symmetry restrictions. The original expression 
of the zero range tensor force formulated by T.H.R. Skyrme (fT4"l) is used (the central and spin-orbit terms are not 
reported); for the sake of convenience, we re- write it here as 



«r(l,2) = 



(<ti • k')fa • k')S(x 1 - x 2 ) + S(xi - x 2 ){(Ti ■ k)fa ■ k) 



U 



T. 
6 



(eri ■ er 2 ) k' 2 S(xi - x 2 ) + 6(xi 



-Ufa ■ k')Sfa - x 2 ){(T 2 -k) - — fa ■ er 2 ) |fc' • 8{x x - x 2 )k \ (1 - P a P T P M ), 



x 2 )k 2 



(Al) 



where k = j=(Vi — V 2 ) and k' — — 57(^1 — V 2 ) respectively act on the right and on the left, er are the spin Pauli 
matrices and P„, P T , Pm are the usual operators which exchange the spin, isospin and spatial coordinates. When 
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computing the Hartree-Fock energy (fTTj) . the (1 — P a P T PM) operator allows, like for the central and spin-orbit terms, 
to account for the exchange without antisymmetrizing the wave- function in the ket. For the considered force, and 
under the hypothesis of no charge mixing, the product of the three exchange operators can be replaced by ±5 qi}q2 
(q = p,n ) for the T and U weighted terms respectively. To emphasize the spatial behaviour of the corresponding 
interaction terms, the tensor parameters arc defined as T — 3t e and U — 3t a in some works. 

Although we deal with a zero range force, it is appropriate to render the dependence on space more explicit in order 
to proceed with the calculation. By representing the a z diagonal basis by the set where lo — 2m s = ±1 

identifies the f, J, eigenvectors, the direct contribution to the Hartree-Fock energy from the first two terms of Eq. (|A1[) 
is computed as 



I, J W 1 / t (Ji ,LJ 2 f W2 



dR 



^c/)*(x' 1 ,uj' 1 )(j)*(x' 2 ,uj' 2 )S(x 1 - x 2 )5(xi ~ cci)5(x 2 - x' 2 ) 



(wj|ffi|wi> • (Vi - V 2 ) (w 2 |tr 2 |w 2 > • (Vi - V 2 ) ^(asi,fa;i)^(x2,W2)}+H.. 



(A2) 



where H.c. denotes the Hermitian conjugate of what preceeds and dR stands for dx' 1 dx' 2 dx\dx 2 (in general, 8{x\ — 
x'i)S(x2 — x' 2 ) ensures the force has a local character). The subscripts i, j summarize the quantum numbers of the 
spinor tpi(x) = Y] ut , 4>i{x,uj)\uj); the isospin is represented through the charge q, among the other quantum numbers. 
In the following, unless differently specified, all the derivative operations apply to the right. 

By inserting the definition of the spin-density (fT2"|) S q (x) = S q (x, x') \ m — m i, where S(x) = ^2 q=p n S q (x) denotes the 
isoscalar density St=o(x), the previous energy contribution can be recast as 



T 
16 



E / dR {[ 2S 

(x 2 , x' 2 )^i,^i, v Sn(xi, x' x ) ~ Vx^S^xx, x' 1 )y 2 ,uS 1/ (x2, x' 2 ) 
-^i, v S^{x l ,x' l )\7 2 ^S v {x 2 ,x 2 ) + H.c. I . 

J J X\—x' 1 —X2—x' 

By inserting the relation (A. 7) of [25j 



and the identity 



]T Sl q H*2, x' 2 ) (Vv,,Vi-,„ + Vi,„Vi,„) S<*> (x 1 ,x' 1 ) = 2S q (x 2 , x' 2 ) ■ G q (x u x[), 



(A3) 



(A4) 



(A5) 



where G(x) = G(x,x')\ x=x > has been firstly introduced in Eq. (fl7|). the energy becomes 

T 



AE? 



J dxidx 2 {2S(xi) • G{x 2 ) - - [Vi • S{x x ) [v 2 ■ S(sb 2 )] + Jbfci) M x *) + 
-\y^^ / i. t ,S u (x 1 )W 2 ^S t ,{x 2 ) + V J IJ , 1/ (x 1 )J 1/IJ ,{x2)\ ■ 

4 £ — ' £ — ' ) xi—xo 



(A6) 



Performing integration by parts twice on the 4th term, after also the S(x\ — x 2 ) function has act, and using, in each 
point of space, the identity 



V - J 2 - - T 2 4- -( tW)< 



(A7) 



for the last term, leads to 



A Ej = T dx 



\s(x) ■ G(x) + 1(V • S{x)) 2 l -J, 2 {x) l -l 2 (x) 



16 



J 2 (x) 



(A8) 



By similarly proceeding for the remaining two terms of the tensor force that are weighted by the parameter T, one 
gets 

AE ?i = ^E /^{[$,(*2,a£)(Vi lM ) 2 S^ 

Z4: — J I L J J xi—x' 1 —X2—x' 2 



2G 



By using the relation (A. 6) of 25 1 

\v 2 + V' 2 )S^(x,x')} = AS^(x)-2T<?\x), 

L ^ J x—x' 

in addition to Eq. (|A4I) , one obtains 



(A9) 



AE 



ii 



T 
24 



J d Xl dx 2 {s{x 2 ) ■ AS( Xl ) - 2S{x 2 ) ■ T( Xl ) - -Va^SUsBajVi^&tei) + 2^ ^(^2)^(^1)} 

}1U 



Once again, due to the zero range nature of the interaction, after one single integration by parts of the last but one 
term, one has 



A Efj = T I dx 



±S(x) ■ AS(x) - is(x) ■ T(x) + if (x) 



the square of the tensor of rank two J is defined, as usual, by J q = ^^(J^u) 2 and, for every point x, the relation 



(A10) 



J q— J-q + Jq + ^i^q* Y 



(All) 



holds. By summing up the two energies AEf and AEf T and taking into account the exchange, the total T-weighted 
contribution to the tensor energy is 



AE° +E 



T 



dx 



\s{x) ■ AS(x) is(x) • T(x) + |(V • S(x) f + S(x) ■ F(x) 



£ 



\s q (x) ■ AS q (x) - ±S g (x) ■ T q (x) + |(V • S 9 (:r)) 2 + S,(x) ■ 



where our equation 



- /" (V • S(a;)) 2 da; = 2 / [flf(x) ■ G(x) + S(sc) ■ F(x)]dx 



(A12) 



(A13) 



has been used. As a matter of fact, the by part integration of the square of the divergence of the spin density reads 
(the integral will be omitted) 

(V-S(x)) 2 = ^2vMx)V v Sv{ic) 

= -^3^X2, x' 2 ) V v ^V 1 ^S u (x 1 ,x' l ) + V 1 . 4J V 1 . u S v (x 1 ,x' 1 ) + 
{xi,x[) + 'V v ^'Vi,uS l/ (xi,x' l )\ 



X\—X\ —X2—X-2 



(A14) 



and one can recognize the structure of the S ■ F and S ■ G terms at the right hand side. 
Concerning the terms weighted by the U coefficient, the first contribution reads 



AE ?n = | £ £ £/ dR^^^x'^^ix'^ 

ij ,UJ\ ,W 2 / ,£^2 jlV 

[^\x 1 ,Lu 1 ) ( f >J (x 2 ,u J2 ) ~ M*l,Ui)*j V \*2,U2j\ }> (A15) 
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where the definition \1/ At (y,a;) = V /t 0(a;, <*})\ m and the short-hand notation 5 — S(x% — X2)5(xi — x' l )8{x 2 — x' 2 ) has 
been introduced. By inserting the definition of the spin density, after a few algebric steps one obtains 



A E? n 



U 



E / dR^S v (x2,x' 2 )Vi' tlJ .Vi ! uS^(xi,x' 1 ) - S7v^S fi (x 1 ,x' 1 )V2 <v S„(x2,x' 2 ) - 

-Vl,vSp(xi, x' 1 )V 2 ',p,S u (x 2 , x' 2 ) + S li (xi,x' 1 )V2,u^2>,(iS u (x2,x' 2 ) \. 

J Xi—X / 1 —X2—X 

By using the relation (|A4|) and the identity 

J2S^(x 2 ,x' 2 ) (Vi^Vi^ + VvViJfiTW^i.xi) = 2S q (x2,x' 2 )-F q ( Xl ,x[), 

fll' 

one is lead to 

A Ef u = - J d Xl dx 2 \2S{x 2 ) ■ F(xi) - - [Vi • S(xi)] [v 2 • S(x 2 )] - Jo (x i)J (x 2 ) 

-^y2Vl. l yS fl (x 1 )V 2 ^S„(x 2 ) - V" J V ^{xx)J^, v {x 2 )\ . 

4 ) Xi— X2 

ill* p,V 

By performing the same steps used to obtain Eq. (|A8[) we get 

U 



(A16) 



(A17) 



Jdx 2S(x) ■ F(x) - i(V ■ S(x)) 2 - J 2 (x) - ]T J^{x)J vlx {x) 



On the basis of Eq. (|A7I) , one finally has 



A E 



in 



2S(x) ■ F{x) -\{V- S(x)) 2 p o 2 (x) + ~J 2 (s) - J 2 (x) 



(A18) 



(A19) 



(A20) 



For the second U contribution of Eq. (|Al|) , the spin structure allows a more compact expression 

U 
~12 



AEfy 



(A21) 



E E E/ ^{*i (m) w,^)^k,-2) 

i,J U^/ ,LJi .W 2 / ,W2 /J, 

o-i • o- 2 <5 (sci, wi)^ (aj 2 , wa) - <f>i (asi, wi)^ rt (» 2 , wa)] } ■ 

This is equal to 

A£f y = -— f dH^{5 1) (s 2l ii!i)Vi t j 1 Vi, (1 S 1/ (ai ) !t' 1 ) - V^lai.xDVj^^y,)},,^^ 

(HI/ 

By using the relations (|A4[) one more time and recognizing the scalar product between the spin and the kinetic density 

E S l q) ( x *> ^Vi' ■ ViSW (xi, X!) = S,(x 2 , x 2 ) • T,(si, si) (A23) 

one obtains 

&Ef v = -— / dx 1 dx 2 {s , (x 2 ) -T(x!) - iE V 2,M^(^2)Vi, M 5^(x!) -E V^i)^^)}^ 

that is 

' ' " ~ (A24) 



AE? V = ~ I ,lx 



S{x) ■ T(x) + ^S(x) ■ AS(x)- J 2 (x) 



By adding together AEfjj and AEf v , one finally gets 



dx T2 



3S(x) ■ F(x) S(x) ■ T(x) - ±S(x)AS(x) - ~(V ■ S) 2 (x) - p 2 (x) + \j 2 (x) - \i 2 {x) 



(A25) 



28 



TABLE IX: Two possible combinations of the D a coupling constants for the direct functional's terms listed in the first row. 
The set (a) corresponds to what is straightforwardly obtained for the direct contributions through the derivation presented 
here, when Eq. (|A13|) is not used in the derivation. The choice (b), on which the result (|A12f) is based, leads to recover the 
expression of [58 1. In both cases, the Galilean invariance of the functional is fulfilled (see Eqs. 



set 


SG 


SF 


(V ■ sf 


(a) 
(b) 


T 

4 




u_ 

\(T+U) 


MT-U) 



The direct plus exchange provides 



AE°+ E 



^Jdx Us(x) ■ F(x) S(x) ■ T(x) \s{x) ■ AS(x) ^(V • S) 2 (x)- 



\j Q 2 (x) + \j*{x)-\l\x) 



3S q (x) ■ F q (x) - S q (x) ■ T q (x) - ~S q (x) ■ AS q (x) - |(V • S q ) 2 (x) 



q—p, 7i 



(A26) 



The sum of the integrand in Eqs. (|A12[) and (|A26|) provides the tensor contribution to the standard Skyrme energy 
density functional, denoted by Uteris in Eq. (fT5|) . 

By defining as D a (E a ) the weight of a generic term a belonging to the direct (exchange) contributions to the energy, 
the relations Z) a =CQ-Cf and D a + E a =C^+Cf hold, where Cy_ 1 are the isoscalar and isovector coupling constants 
employed in some other works. By rearranging the result in terms of the isoscalar and isovector densities, the sum of 
the expressions (49a) + (49b) of [58[ is recovered. Otherwise, when the S ■ G terms from (IA8[) are retained, the S ■ F 
and (V • S) 2 coupling constants in (|A12I) change according to the combination (a) of Tab. IIX[ in such a way, one 
obtains, in the proton- neutron formalism, the more general formulation (|15p . where A a = D a + E a and B a = D a . 
Since E a = ^fD a for the T (-) and U (+) contribution respectively, only the U-weigthed terms contribute to A a . 



Appendix B: Densities and currents of the 
Skyrme-tensor EDF 

In this Appendix, the expressions of the densities 
and currents entering the Skyrme energy density func- 
tional ([L3)) are provided. For the sake of semplicity, 
we take implicit the sum over the wave-functions index, 
which includes all the quantum numbers with the excep- 
tion of the spin projection oj, which is separately written. 
When boosting the Hartree-Fock single-particles, they 
gain an imaginary part. We adopt the definitions of 
Tab. [Xl where 0^ = 4>(x, ui), k = x,y,z and u> — 2m s 
■ 1 . and of Tab. IXI1 where the second order terms are 
listed. In the following, all the relations hold if consid- 
ering isoscalar, isovector, neutron or proton densities. 

The densities which are even under the time-reversal 
operation Tk ■ <pi(x,oj) —> —u>(f)*(x, —lu) (Tk = —icjyKo, 
where Kq denotes the complex-conjugation), include 

• the particle density 



TABLE X: Conventions (I) adopted in the current Appendix. 



4>S 












n 




n 






V(f> u 






k 













• the kinetic density 

t(x) = {V-V'p(x,x')} x = x ,; (B2) 

• the spin current pseudo-tensor J 

J (x) = j. [(V - V) ® S(x, x')] x=x , . (B3) 
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TABLE XI: Conventions (II) adopted in the current Ap- 
pendix. 



The 9 components can be implemented as 

j kx = E 

UJ.LO'—^zl^UJ^Oj' 

L=R,I 



+3 [<a;*+-c* s 



(B4) 
(B5) 



(B6) 

The time-odd densities are 
• the spin density 

S(x) = E P(x',u')ct>(x,tj)(uj'\*\cj)\ x=x , , (B7) 



where each component can be worked out as 
S x = 2[<pK<P* + cj> I + cj> I _}=2$l[<j>* + <P_] 
S y = 2[<f> R _<f> I _ - </>*(/>'+] = 23[«^_] 

Sz =52 i<i>+<t>+ - 4> L -i> L -\ = m>*+4>+ - <? 

L=R,I 

• the momentum density 



(B8) 
(B9) 
-IB10) 



j(x) = --{(V-V')p(x,x')} a 
with components equal to 

h - 52 [ W 

• the spin kinetic density 

T{x) = [V-V'S{x,x')] x=x ,, 

where 

T x = 2[^l^ R + ^^l] 
= 23?[*_#+] 

23?[*_*+] 

52 



(Bll) 



T — 



L=R.I 



(B12) 
(B13) 

(B14) 
(B15) 

(B16) 



the pseudo-vector F(a;) density 



F(x) - - {[(V ® V) + (V ® V')] 5(x, *')}„=„' . 



It can be conveniently represented by splitting each 
fc=x,y,z component in three contributions 



F k (x)= 52 Fi k '\x), 

k'—x,y,z 

where the single terms read 
= 2R[*+*"] 

= - *+] 

E [*J' + *x , - + *J'-*x' + ] 

2KE' / '~vJ/ fl ' + - ii/ fl '-\l/ i '>+l 

i v v y v J 

-2S[*;*+] 

E [*j >+ *^ + - 



- 1 z 



= fr( x ) I 

= E K*J ,+ ) a - (^'") 2 

L=R,I 



(B17) 



(B18) 



(B19) 



(B20) 



(B21) 

(B22) 
(B23) 

(B24) 

(B25) 
(B26) 

(B27) 



• the pseudo- vector C?(aj) density 



= 2 {[(V ® V) + (v ® v)] s(sb, *')}„=„/ ; 
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by proceeding similarly to F(x), we obtain 

G ( x x) = £ [ft A x <p L _ + J> L _ A x 4> L + ] 

L=R,I 

= A, </>_ + <j)*_ A x 0+] (B28) 

G x v) = o'M'!.; ( ,/«!'»• o' >I'»-' 

= -0**+] (B29) 

= £ [<ft*Zi + -4>^i-] 

L=R,I 

= -01*-] (B30) 



4 X) = E [0+*£r + 0-^ 



= (B31) 

G(») - 0* A y 7 _ - ^ A, ^' + 
0i A y ^ - K A„ 
= A y 0_ - 0*_ A,, 0+] (B32) 

G ^ = (B33) 



= G^|^ 2 (B34) 
= G x y\^ z (B35) 

g{ z) = E A * ^+ - A * 

L=R,I 

= $t[<t>* + A z cf> + - <{>*_ A z <j>_\. (B36) 



